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ABSTRACT 


This  thesis  examines  the  system  implementation  considerations  of  a  resonance- 
based  radar  target  classification  system.  The  basis  of  the  system  is  the  aspect  and 
excitation  independent  property  of  electromagnetic  scattering  from  a  conducting 
body.  Such  a  system  consists  of  two  components:  pole  extraction  and  annihilation 
filtering.  The  algorithms  investigated  here  for  these  purposes  are  the  Cadzow- 
Solomon  pole  extraction  algorithm  and  the  K-Pulse  annihilation  filter.  Additionally, 
an  aspect-dependent  annihilation  filter  based  on  an  inverse  ARMA  model  is 
introduced.  The  procedures  are  applied  to  noise  polluted  synthetic  data,  as  well  as 
scattering  data  collected  for  a  thin-wire  and  silver  coated  1/72  scale  model  aircraft. 
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1.  INTRODUGTIGN 

Modern  antiair  warfare  has,  with  the  advent  of  radar  in  World  War  II,  become 
a  process  of  detecting  aircraft  and;  identifying  them  as  friend  dr  foe  at  sufficient 
ranges  to  allow  engagement  with  appropriate  weapon  systems.  Radar  developriient 
has  historically  emphasized  the  first  element  of  this  process:  detection.  As  modern 
air  warfare  has  become  more  complex,  dynamic,  and  lethal,  identification  has  become 
as  critical  as  detection. 

In  the  past,  identification  has  been  facilitated  by  active  measures  on  the  part  of 
friendly  aircraft  to  distinguish  them  from  enemy,  i.e.  IFF  transponders,  specific  flight 
routes,  etc.  AJl  these  are  subject  to  exploitation  by  an  adversary.  This  has  given  rise 
to  development  of  non-cooperative  target  recognition  (NCTR),  the  ability  to  identify 
an  enemy  aircraft  despite  its  efforts  to  deceive. 

One  approach  to  NCTR,  as  discussed  in  this  thesis,  is  identification  based  on 
the  natural  complex  resonances  of  the  aircraft  as  a  scatterer  of  incident 
electromagnetic  (radar)  energy.  In  1971  Baum  showed  that  an  incident 
electromagnetic  wave  creates  a  current  response  on  a  body  with  certain  complex 
natural  resonances  dependent  on  its  physical  configuration  [Ref.  1].  Moffat 
and  Mains  proposed  that  a  target  scatters  electromagnetic  energy  with  resonances 
that  are  independent  of  incident  polarization  or  angle  [Ref.  2].  Morgan 
confirmed  this  through  a  rigorous  theoretical  description  of  induced  currents  and  the 
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resultant  scattered  field;  [Ref.  3].  The  problem,  then,  is  to  determine  these 
complex  poles  for  a  target  of  interest  and  subsequently  use  them;  to  identify  this 
target. 

The  first  step  in  using  complex  natural  resonances  for  target  identification  is 
extraction  of  these  poles  from  scattered  responses  of  actual  targets  (either  modeled 
or  full  size)  [Ref.  4].  These  poles  are  then  used  to  create  annihilation  filters 
for  each  potentially  observed  target  [Ref.  5].  A  given  filter  will,  in  principle, 
produce  a  near  null  late-time  response  to  its  matched  target’s  scattered  signal, 
independent  of  aspect  or  polarization. 

The  objective  of  this  thesis  is  to.  examine  these  processes  in  order  to  define  the 
overall  system  characteristics  and  operational  considerations  of  a  radar  target 
identification  system  based  target  natural  resonsnces.  This  analysis  will  use  the 
application  of  signal  processing  algorithms  to  both  synthetic  and  actual 
electromagnetic  scattering  data.  This  thesis  has  been  conducted  in  conjunction  with 
research  in  pole  extraction  by  Murphy  [Ref,  6],  whereby  the  data  generated 
in  that  work  was  incorporated  in  the  filters  discussed  in  this  effort. 

A.  BACKGROUND 

The  ability  to  extract  poles  and  use  them  to  identify  a  target  is  dependent  on 
the  nature  of  electromagnetic  scattering  from  a  conducting  body.  Consider  a 
generalized  electromagnetic  pulse  with  electric  and  magnetic  field  components 
and  B{r,t)  respectively,  travelling  in  a  direction:  p  and  incident  upon  a  finite-sized. 
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perfectly  conducting  body  iniree  space.  Currents  are  induced  on  the  body  which  are  . 
described  by  the  magnetic  field  integral  equation  [Ref.  7], 

Jir,t)  =  2nxH\r,t)  +  ff  K(ry,t)‘J(r',t-^^^)dS  (1) 

c 

where  «  is  the  outward  unit  normal  vector,  /  is  the  surface  current  density,  and  K  is 
a  dyadic  Green’s  function  kernel.  Here  the  surface  integral  is  a  principle-value  type, 
excluding  the  point  r  =  r' .  The  cross  product  term  represents  the  physical  optics 
contribution  to  the  induced  currents  and  the  surface  integral  represents  the  feedback 
of.currents  on  each  point  on  the  object  due  to  the  currents  at  all  other  points.  This 
feedback  mechanism  forms  the  physical  basis  for  natural  resonances. 

After  the  incident  pulse  leaves  the  target  body,  JFT  in  Equation  (1)  goes  to  zero; 
the  source-free  current  modes  that  remain  "ring"  at  the  natural  resonance  frequencies 
of  the  object.  These  currents  are  described  by  the  expression /„(r)exp(5,/);  where 
=  (7„  +  are  the  frequencies  which  are  dependent  only  on  the  physical 
properties  of  the  scatterer.  These  complex  frequencies  are  described  in  the  Laplace 
domain  by  their  damping  constant,  a,  in  Nepers/sec,  and  their  frequency,  to,  in 
radians/sec.  An  infinite  number  of  these  resonant  modes  are  inherent  in  any 
candidate  scattering  object,  but  only  those  in  the  bandwidth  of  the  incident  pulse  are 
significantly  excited.  [Ref.  7] 
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The  induced  currents  on  the  scattering  body,  J,  give  rise  to  a  far-field 
distribution  in  the  backscatter  direction,  -p,  (the  direction  o£  interest  for  monostatic 
radar  applications)  described  by 

H^(-rp,t)  =  ^  (2) 

47:cr  c 


This  equation  is  applicable  during  the  entire  time  the  incident  pulse  is  traversing  the 
target  in  the  direction  p  as  well  as  after  it  has  gone  beyond  it.  When- the  specific 
current  distribution  described  in  Equation  (1)  is  substituted  here,  the  result  is 


H^(-rp,f)  -  u(t-r/c) 


J  H^(rrp,t)  + 


Rt‘0 


(3) 


where  Hp^  is  a  physical  optics  term  dependent  on  the  incident  direction,  p,  relative 
to  the  object  shape,  and  the  summation  term  is  due  to  the  source-free  currents 
behind  the  advancing  wavefront.  The  field  due  to  these  currents  resonates  at  the 
same  frequencies  as  the  currents  themselves.  However,  during  the  time  the  pulse  still 
resides  somewhere  on  'the  target,  the  coefficients  H„  are  varying.  After  the  pulse 
leaves,  the  aspect  dependent  physical  optics  term  vanishes  and  the  coefficients  of  the 
scattered  resonant  frequency  modes  become  constant.  This  instant  defines  the 
transition  from  "early-time"  to  "late-time"  in  the  scattered  field.  For  backscattering 
this  transition  occurs  after  a  time  given  by  T„  =  2D/c  +  t,  as  measured  from  when 
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the  pulse  first  impinges  on  the  target,  where  D  L  the  length  of  the  object  in  the 
propagation  direction,  c,  the  speed  of  light;  and  t,  the  incident  pulse  width.  [Ref. 

7] 

B.  PROBLEM 

Application  of  these  characteri.stics  to  target  identification  involves  tv/o  separate 
steps:  pole  extraction  and  annihilation  filtering.  The  transfer  function  of  an 
electromagnetic  target  contains  a  complicated  set  of  poles  that  can  be  analytically 
predicted  only  for  the  simplest  objects,  e.g.,  a  thin-wire  or  a  sphere  [Ref.  8]. 
Therefore,  it  will  be  necessaiy  to  extract  the  poles  from  complex  targets,  such  as 
aircraft,  through  signal  processing  of  their  scattered  response  [Ref.4].  A  number  of 
algorithms  have  served  as  candidates  for  this  process  [Ref.  4],  including  the 
Kumaresan-Tufts  and  Cadzovz-Solomon  algorithms  compared  by  Larison  in 
[Ref.  9].  These  will  be  discussed  further  in  Chapter  III. 

After  the  poles  are  determined  for  the  targets  of  interest,  they  will  be 
incorporated  in  a  bank  of  annihilation  filters;  each  filter  corresponding  to  a  target. 
Scattering  responses  from  unknown  targets  (radar  returns  in  actual  application)  are 
passed  to  all  the  filters.  The  annihilation  filter  with  the  least  output  should  correlate 
the  unknown  target  to  a  previously  cataloged  response,  thus  identifying  that  target. 
This  is  the  basis  of  the  Resonance  Annihilation  Filter  (RAF)  concept  discussed  in 
[Ref.  7]  and  [Ref.  10].  Further  analysis  of  this  process  and  its  implications 
for  actual  implementation  are  discussed  in  Chapter  IV. 
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G.  HISTORY 


As  stated,  the  initial  work  in  the  areaiof  natural  resonance  response  was  derived 
from  Baum’s  investigation  of  currents  generated  as  airesult  of  electromagnetic  pulse 
(BMP)  [Ref.  1].  Tesche  developed  the  concept  further  into  the  specifics  of  the 
Singularity  Expansion  Method;(SEM)  applied  to  generating  the  time  response  of  a 
thin  wire  scatterer  analytically  by  examining  its  poles  and  residues  [Ref.  11^. 
Morgan  points  out  in  [Ref.  3]  and  [Ref.  7]  that  the  SEM  description  of  the  scattered 
signal  is  applicable  only  in  late-time,  the  time  after  the  incident  pulse  leaves  the 
object.  As  a  result  of  this  limitation,  earlier  schemes  for  pole  extraction  based  on 
Prony’s  method  [Ref.  12]  were  limited  by  noise  and  clutter  [Ref.  13J. 

Since  that  time,  effort  has  been  directed  by  many  researchers  into  improved 
pole  extraction  algorithms.  Chong  incorporat  cd  ihe  early-time  in  using  a  direct  non¬ 
linear  least-squares  minimization  [Ref.  14],  but  this  approach  is 
com.putationally  intensive  and  requires  much  a  priori  knowledge  of  the  signal  [Ref. 
4].  The  Kumaresan-Tufts  method  of  pole  extraction  is  limited  in  applicability  to  the 
late-time  only,  but  relieves  some  of  the  burden  of  a  priori  knowledge  of  the  natme 
of  the  scattered  re.sponse  [Ref  15]  and  improves  the  noise  handling 
capability  of  Prony’s  method.  The  Cadzow-Soior.  ..  n  technique,  a  different  e;dension 
of  Prony’s  method,  makes  use  of  both  the  e?;  ly  and  late-time,  allowing  it  to  use  as 
much  of  the  available  eneigy  as  desired  in  the  return  signal.  This  algorithm, 
introduced  in  [Ref.  16]  and  applied  to  pole  extraction  in  [Ref.  4],  [Ref.  9] 
and  [Ref.  6],  does  require  input  of  the  transmitted  signal  along  with  the  return 
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scattering  signal,  but  holds  the  most  promise  of  the  considered  extraction  methods 
in  the  presence  of  noise. 

The  annihilation  filtering  concept,  described  above,  has  been  attempted 
primarily  with  two  basic  algorithms.  I  i-  ulse  algorithm,  based  oma  difference 
equation  representation  of  the  filter,  cn  developed  by  Morgan  [Ref.  7]  and 
Dunavin  [Ref.  17]  from  the  orip.i'ial  v/ork  of  Keniiaugh  [Ref.  18). 
The  E-Puise  technique,  developed  by  [Ref.  19],  uses  integral  operators 

to  cancel  the  resonances  in  late-time.  Simon  compared  these  algorithms  in  [Ref.  10: 
and  found  similar  results.  This  thesis  will  further  investigate  the  K-Pulse  annihilation 
filter,  and  will  extend  it  to  a  recursive  version  which  works  on  early-time  energy  as 
well. 
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II.  EXPERIMENTAL  SIGNALS 

Both  synthetic  and  actual  signals  were  used  to  analyze-ths  performance  oLpole 
extraction  and  annihilation  filtering.  Two  types  of  syntheticrsignals  were  used.  First, 
synthetic:  signals  with  predetermined  complex  poles  were  used  to  valriate  extraction 
anc  annihilation  software  rnd  to  gain  insight  into  the  mechanisms  by  which  these 
algorithms  worked.  Secondly,  poles  were  extracted  from  syht  utically  generated  thin- 
wire  scattering  data  generated  by  Morgan’s  thin-wire  integral  eq.idtion  program 
[Ref.  20],  and  used  in  annihilation  filters.  For  this  set  of  data,  the  poles 
were  not  known  in  advance.  After  the  software  was  tested  on  synthetic  data,  it  was 
applied  to  actual  scattering  data,  gathered  by  Walsh  [Ref.  21].  This  data  was 
measured  for  a  thin  wire  and  several  scale  model  aircraft  using  the  NFS  Transient 
Electromagnetic  Scattering  Luooratory  (TESL).  In  the  development,  of  the  first  set 
of  synthetic  signals,  as  many  parameters  as  possible  were  selected  to  iiMtch  those 
anticipated  in  the  TESL  measured  data. 

A.  ARMA  SYNTHETIC  SIGNALS 
1.  Scatterer  Model 

As  discussed  in  the  previous  chapter^  the  scattered  field  consists  of  a 
physical  optics  term  due  to  the  driving  incident  wave  and  a  feedback  portion 
representing  the  complex  natural  resonances.  These  characteristics  suggest  an 
autoregressive  moving  average  (ARMA)  or  recursive  structure  for  modeling  this 
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signal.  Such  a  structure  was  realized  with  a  recursive  difference  equation  signal 
generator, 


yin) 


k^N  k=L 


=  E  ^tyi^-k)  +  '^b^in-k) 

i=l  i.O 


(4) 


where  a:(«)  is  the  digital  version  of  the  incident  pulse,  >»(n)  is  the  scatterer  response 
and  and  are  coefficients  determined  by  the  complex  natural  resonances  and 
initial  conditions,  respectively  [Ref.  22:pp.  81-84].  After  selection  of  a 
sampling  frequency  for  the  synt’ietic  signals,  the  coefficients  were  determined  through 
conversion  of  a  simple  set  of .  -plane  poles  to  the  complex  z-plane,  followed  by 
multiplication  of  rational  polynomial  terms  representing  each  pole  pair  to  arrive  at 
an  overall  ratio  of  polynomials:  the  z-trahsform  verf.ion  of  the  scatterer’s  transfer 
function.  The  coefficients  from  the  numerator  and  denominator  appear  in  Equation 
(4)  as  bf.  and  a^.,  respectively.  Additional  inputs  were  a  complex  coefficient  of  the  z- 
plane  poles  to  allow  scaling  or  phase  differences  betv/een  poles.  Program  listings  for 
the  coefficient  generator  and  the  recursive  signal  generator  appear  in  Appendices  A 
and  B. 

The  poles  were  selected  in  the  ^-plane  to  coincide  generally  in  bandwidth 
with  those  obtainable  by  actual  scattering  measurements.  As  such,  three  sets  of  ten 
pole  pairs  each  were  selected,  lying  at  integer  values  between  1-10  GHz.  The  three 
set:  represent  different  levels  of  damping  (Q  factor)  and  are  related  in  that  they  all 


have  the  same  modulus  when  plotted  as  vectors  in  the  5-plane.  For  each  set  of  poles, 
the  frequency  and  damping  coefficient,  and  cr^,  are  related  by  a  daniping  factor 


(5) 


where 

kj  =  0.5  (low  Q/highly  damped) 

ko  =  0.7  (medium  Q/moderately  damped) 

k^  =  0.8  (high  Q/lighOy  damped)  [Ref.  10]. 

These  poles  are  plotted  in  the  5  and  z-planes  respectively  in  Figures  l=and'2  and  are 
listed  in  Tables  1,  2  and  3.  The  units  of  and  cd^  are  GHz,  GNp/sec  and 
GRad/sec,  respectively,  The  resulting  ARMA  coefficients,  Of.  and  bf.,  art  in  Table  4. 


TABLE  1.  LOW  Q  SYNTHETIC  POLES 


POLES 

"p 

RESIDUES 

1 

-0.6892 

6.2474 

1 

0 

2 

-13784 

12.4948 

1 

0 

3 

-Z0676 

18.7422 

1 

0 

4 

-Z7568 

24.9895 

1 

0 

5 

-3.4460 

31.2369 

1- 

0 

6 

-4.1352 

37.4843 

1 

0 

7 

-4.8244 

43.7317 

1 

0 

8 

-5.5136 

49.9791 

1 

0 

9 

-6.2028 

56.2264 

1 

0 

10 

-63919 

62.4738 

1 

0 

10 


Im  s  (Gigarad/sec) 


TABLE  2.  MEDIUM  Q  SYNTHETIC  POLES 

POLES  RESIDUES 


«/> 

1 

-03562 

6.2752 

1 

0 

2 

-0.7124 

123504 

1 

0 

3 

-1.0687 

183256 

1 

0 

4 

-1.4249 

25.1007 

1 

0 

5 

-1.7811 

313759 

1 

0 

6 

-Z1373 

37.6511 

1 

0 

7 

-Z4935 

43.9263 

1 

0 

8 

-23498 

50.2015 

1 

0 

9 

-3.2060 

56.4767 

1 

0 

10 

-33622 

6Z7518 

1 

0 

fp 

TABLE  3.  HIGH  Q  SYNTHETIC  POLES 

POLES  RESIDUES 

4. 

1 

-0.1625 

6.2832  1 

0 

2 

-03250 

123664  I 

0 

3 

-0.4876 

183196  1 

0 

4 

-0.6501 

25.1327  1 

0 

5 

-03126 

31.4159  1 

0 

6 

-0.9751 

37.6991  1 

0 

7 

-1.1376 

43.9823  1 

0 

8 

-13002 

50.2655  1 

0 

9 

-1.4626 

563437  1 

0 

10 

-1.6252 

623319  I 

0 
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The  sampling  frequency  of  25.6  GHz,  512  samples  in  20  nanoseconds,  was 
used  in  the  creation  of  the  synthetic  signals.  This  is  the  minimum  "power-of-two" 


TABLE  4.  AJRIVU  COEFFICIENTS  FOR  SYNTHETIC  SCATTERING 

Low  Q  Medium  Q  High  Q 


k 

•k 

h 

•k 

h 

•k 

h 

0 

1.0 

20.0 

1.0 

20.0 

1.0 

20.0 

1 

-3.77 

-71.59 

-3.54 

-67.32 

-3.46 

•65.65 

2 

8.47 

152.61 

8.08 

145.36 

7.96 

143.26 

3 

-14.27 

-242.52 

-14.09 

-239.57 

-14.14 

-240.43 

4 

20.02 

320.34 

21.00 

335.97 

21.66 

346.64 

5 

-24.51 

-367.64 

-27.61 

414.21 

-29.40 

441.01 

6 

nsn 

379.04 

33.15 

464.08 

36.58 

512.05 

7 

■27.40 

-356.18 

-36.69 

476.94 

42.03 

-546.42 

8 

25.80 

309.55 

38.02 

456.23 

45.32 

543.92 

9 

-22.70 

•249.71 

-36.96 

406.56 

45.90 

-504.90 

10 

18.84 

188.36 

34.00 

339.98 

44.02 

440.24 

11 

-14.72 

-132.44 

-29.50 

•265.52 

•39.85 

-358.65 

12 

10.89 

87.13 

24.28 

194.28 

34.22 

273.76 

13 

-7.58 

-53.04 

-18.80 

■131.63 

-27.64 

-193;49 

14 

4.99 

29.92 

13.75 

82.51 

21.10 

126.40 

15 

-3.04 

-15.22 

-9.33 

46.66 

-14.89 

-74.45 

16 

1.74 

6.98 

5.90 

23.60 

9.77 

39.09 

17 

-0.89 

-2.67 

-3.33 

-9.99 

-5.72 

-17.17 

18 

0.42 

0.84 

1.69 

3.39 

2.99 

5.98 

19 

-0.15 

-0.15 

-0.67 

■0.67 

-1.21 

-1.22 

20 

0.05 

0.00 

0.22 

0.00 

0.38 

0.00 

number  of  samples  allowable  in  a  20  nanosecond  collect  time  by  the  Nyquist 
requirement  applied  to  the  frequencies  of  interest,  one  to  ten  GHz.  The  requirement 
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to  limit  the  sampling  bandwidth  as  much  as  possible  is  derived  from  the  use  of  the 
annihilation  filter,  discusseo  Chapter  IV. 

2.  Incident  Pulse 

After  determination  of  the  coefficients  of  the  synthetic  scatterer  transfer 
function,  these  were  used  in  the  difference  equation,  Equation  (4),  above.  The  input 
waveform,  x(n),  was  chosen  to  be  a  double  Gaussian,  in  order  to  match  that  used  in 
the  actual  scatterer  data  collection.  The  double  Gaussian  waveshape  consists  of  a 
broad,  slow  Gaussian  pulse  subtracted  from  a  narrow,  rapid  Gaussian  pulse.  The 
narrow  pulse  determines  the  high  frequency  roll-off;  the  broad  pulse  determines  the 
low  frequency  roll-off.  The  incident  pulse  used  here  has  ten  percent  narrow  and  wide 
pulse  widths  of  0.15  and  0.3  nanoseconds,  respectively,  and  is  pictured  in  Figure  3. 
The  frequency  domain  representation  is  in  Figure  4.  The  final  results  of  the  double 
Gaussian  input  to  the  ARMA  model  is  shown  for  all  three  pole  sets  (high  Q,  medium 
Q,  and  low  Q)  in  Figures  5,  6,  and  7  and  Figure  8  contains  the  spectrum  of  the 
medium  Q  target. 

3.  Noise  Modeling 

The  synthetic  signals  were  polluted  by  simulated  Gaussian  noise  at  various 
signal-to-noise  ratio  levels:  20  dB,  15  dB,  10  dB  and  7  dB.  The  noise  was  introduced 
by  a  noise  pollution  routine,  Appendbc  C,  which  uses  a  random  number  generation 
algorithm  based  on  subroutines  of  the  Fortran  Scientific  Library.  This  routine 
generates  random  numbers  under  a  uniform  distribution.  The  Central  Limit 
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Magnitude 


.  Frequency  (GHz) 

Figure  4.  Frequency  spectrum  of  double  Gaussian 
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Figure  5.  High  Q  synthetic  scattering 
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Figure  6.  Medium  Q  synthetic  scattering 
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Figure  7.  Low  Q  synthetic  scattering 
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Magnitude 


Figure  8.  Frequency'  spectrum  of  medium 
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Theorem  was  then  invoked  by  calling  this  routine  1000  times  for  each  additiye  noiie 
point  to  ensure  the  resultant  noise  was  Gaussian  distributed.  [Ref.  6];  Figures  9  and 
10  show  the  high  Q  waveform  of  Figure  5  with  the;addition  of  noise  of  20  dB  and  7 
dB,  respectively. 

B.  THIN-WIRE  INTEGRAL  EQUATION  SIGNALS 

The  concept  of  extracting  nnfoiown  poles  and  then  incorporating  them  into 
filters  was  validated  on  the  signals  generated  by  the  Time-Domain  Integral  Equation 
(TDiE)  Thin-Wire  Scattering  Program,  due  to  Morgan  [Ref.  20].  This  program  uses 
a  record  length  of  960  samples  in  a  20  nanosecond  collect  time.  The  scattering 
simulated  excitation  by  the  same  double  Gaussi«>n  pulse  with  approximately  a  1-10 
GHz  bandwidth.  The  program  generates  wire  scattering  data  for  any  combination 
of  incident  and  reflection  angles.  Backscatter  data,  generated  at.  angles  of  175, 150, 
110,  and  90  degrees  measured  relative  to  the  long  axis  of  the  wire,  were  used  to 
effect  a  monostatic  situation.  This  data  was  similarly  noise  polluted  at  the  same 
signal-to-noise  ratios.  A  representative  case  for  the  noiseless  TDIE  thin-wire  data  is 
in  Figure  11  and  its  spectrum  is  in  Figure  12. 

C.  LABORATORY  SIGNALS 

Actual  scattering  signals  used  for  pole  extraction  and  annihilation  filtering  were 
obtained  from  the  NPS  Transient  Electromagnetic  Scattering  Laboratory.  As  stated, 
the  parameters  determined  by  the  instrumentation  of  the  laboratoiy  were  used 
throughout  the  experiments  for  both  synthetic  and  actual  signals  as  much  as  possible. 
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Frequency  (GHz) 

Figure  12.  Spectrum  of  TDIE  thin-wire  307150®  datn 
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This  included  the  scattering  collect  time  of  20  nanoseconds  and  the  excitation 
bandwidth  of  1-10  GHz  due  to  the  incident  pulse?  double  Gaussian  waveshape.  The 
sampling  frequency  serves  as  the  notable  exception  to  this  matching.  The  TESL 
samples  1024  times  in  the  20  nanosecond  collect  time  for  a  sampling  frequency  of 
51.2  GHz.  This  twice  oversampling,  while  allowing  more  flexibility  in  the  use  of  the 
TESL,  acted  as  a  limiting  factor  for  subsequent  annihilation  filtering.  As  a  result,  the 
collected  scattering  data  from  the  lab  was  collapsed  to  the  25.6  GHz  sampling 
fi’equency,  mentioned  above,  by  averaging  each  pair  of  points  in  the  record  before 
annihilation  filtering.  The  TESL  provides  a  nominal  signal-to-noise  ratio  of 
approximately  26  dB  [Ref.  21]. 

1.  Deconvolution  Postprocessing 

The  double  Gaussian  waveshape,  as  described  above,  is  not  actually 
transmitted  in  the  TESL.  Instead,  a  noise-like  but  highly  repetitive  signal,  with 
approximately  the  same  bandwidth  is  transmitted  and  the  response  of  the  target  to 
the  specific  double  Gaussian  waveshape  is  estimated  through  deconvolution  post¬ 
processing.  This  allows  the  target  to  be  subjected  to  an  analytically  precise  signal  like 
the  double  Gaussian  without  requiring  a  tremendously  sophisticated  transmitter 
apparatus.  Additionally,  the  free-space  response  of  the  target  is  more  nearly 
obtained  by  removing  the  effects  of  the  chamber  and  antennas.  The  details  of  this 
process  are  in  [Ref.  21].  Targets  used  in  the  scattering  ex'periments  included  a  0.1 
meter  thin-wire  target  and  silver  coated  1/72  scale  model  aircraft. 
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2.  Measured  Thin  Wire  Target 

The  transition  to  real-world  scattering  data  was  first  made  with  scattering 
data  collected  in  the  TESL  on  a  0.1  meter  long,  2.36  millimeter  diameter  thin  metal 
wire  [Ref.  21].  This  wire  data  was  collected  monostatically  at  30  degrees,  45  degrees 
and  broadside  to  the  wire.  The  time  and  frequency  domain  data  for  the  30  degree 
case  are  in  Figures  13  and  14.  Notice  that  the  FFT  of  this  target’s  scattering  is  much 
more  complicated  than  for  the  synthetic  TDIE  thin  wire.  The  spectrum  in  Figure  14 
displays  a  significant  resonance  at  approximately  3  GHz,  with  less  excited  poles  at  1.5 
and  4  GHz.  As  the  apparent  length  of  the  wire,  D,  at  this  angle  is  0.05  meters,  so 
these  resonances  correspond  roughly  to  modes  at  »  2D,  3D,  and  1.5D,  respectively. 
This  demonstrates  the  complicated  interaction  between  the  excitation  and  resonances 
of  the  target. 

3.  Model  Aircraft  Targets 

The  bandwidth  of  the  TESL,  1-12  GHz,  is  matched  to  the  scaling  factor 
of  the  model  aircraft  targets.  Table  5  contains  the  significant  dimensions,  both  actual 
and  scaled,  of  the  aircraft  used  to  collect  scattering  data.  Comparison  of  the  scale 
model  lengths  to  the  wavelengths  of  the  TESL  excitation,  2.5  to  30  cm,  indicate  that 
large  fractions  of  wavelengths  can  be  contained  in  the  dimensions  of  the  aircraft 
models.  This  allows  for  the  lowest  order,  i.e.  the  strongest  available,  modes  to  be 
excited  in  the  model  structure. 
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Time  (nanoseconds) 


Figure  13.  Measured  thin-wire  scattering  at  307150° 
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TABLE  5.  DIMENSIONS  OF  AIRCRAFT  TARGETS  [Ref.  23,  24] 


Target  number 

1 

2 

3 

4 

5 

Overall  length 

actual  (m) 

19.43 

15.03 

17.07 

12.20 

15.88 

model  (cm) 

26.97 

20.88 

23.71 

16.94 

22.06 

Wingspan 

actual  (m) 

13.05 

9.45 

11.43 

10.96 

13.95 

mode!  (cm) 

18.13 

13.13 

15.88 

15.22 

19.38 

The  aircraft  data  was  collected  by  Walsh  at  six  different  aspects:  nose^on, 
tail-on,  broadside,  30  degrees  off  the  nose,  and  from  above  with  both  the  fuselage 
and  wingspan  parallel  to  the  incident  electric  field  [Ref.  21].  Target  No,  2  scattering 
data  at  thirty  degrees  aspect  angle  is  in  Figure  15  with  its  spectrum,  Figure  16. 

D.  SUMMARY 

Each  of  these  target  sets;  the  ARMA  synthetic  scatterer,  the  TDIE  thin-wire 
scatterer,  the  actual  thin  wire,  and  the  model  aircraft;  each  represent  an  increasingly 
more  complicated  and  challenging  case  for  the  pole  extraction  and  annihilation 
filtering  algorithms.  This  is  especially  evident  when  reviewing  their  various  frequency 
response  plots. 
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Figure  15.  Scattering  from  aircraft  target  No.  2,  30®  aspect 
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III.  POLE  EXTRACTION 

The  pole  extraction  algorithm  considered  in  this  thesis  is  due  to  Cadzow  and 
Solomon  [Ref.  16].  Other  algorithms  have  been  used  and  are  briefly  discussed  in 
Chapter  I.  More  detailed  analysis  of  the  other  algorithms  can  be  found  in  [Ref.  4]. 

A.  CADZOW-SOLOMON  ALGORITHM 

The  scattered  signal  structure  discussed  in  Chapter  I  and  displayed  in  Chapter 
II  consists  of  an  early-time;  and  late-time  separated  by  the  point  where  the  incident 
signal  ceases  to  excite  the  scatterer.  In  early-time,  the  resonances  are  just  being 
initialized,  so  their  amplitudes  are  small  compared  to  the  driver.  However,  the  end 
of  early-time  marks  the  point  where  they  have  their  largest  instantaneous  power.  The 
Cadzow-Solomon  algorithm  takes  advantage  of  this  fact  by  attempting  pole  extraction 
from  a  data  record  containing  both  early  and  late-time  data  segments. 

1.  Applicability 

The  Cadzow-Solomon  method  seeks  to  identify  the  unknown  coefficients 
and  for  an  autoregressive  moving  average  (ARMA)  process  given  by  Equation 
(4),  which  is  repeated  here, 

yin)  =  E  nf^in-k)  +  ]E  h^in-k\  (<») 

i=l  t=o 
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by  analyzing  the  system  input,  :c(/i)  and  output,  _y(n)  [Ref.  16],  Identificatiomof  these 
coefficients  will  allow  for  extraction  of  the  poles  and  zeros  of  the  transfer  function 
which  has  been  acted  on  byx(n)  to  produce  y(n),  specifically, 


H(z)  = 


^0 

1  +  a^z'^  + 


•  •  + 


(7) 


2.  Methodology 

The  algorithm  does  not  require  specific  knowledge  of  the  system  order,  N 
and  L  in  Equation  (7),  as  it  uses  an  "over-estimation"  of  this  order  to  account  for 
noise  and  extract  the  poles.  Norton  applied  a  technique  used  in  the  Kumaresan-Tufts 
algorithm,  whereby  the  relationship  in  Equation  (4)  is  written  in  reverse  time  order. 
This  forces  the  desired  poles  to  be  extracted  outside  the  unit  circle  on  the  z-plane, 
while  the  extraneous  poles,  introduced  by  "over-estimation"  to  allow  for  noise,  are 
constrained  to  lie  inside  the  unit  circle  [Ref.  4].  This  reverse  time,  non-causal  form 
of  Equation  (4)  is 

k=N  k=L 

y(n)  =  5^ajty(N+«-A:+l)  +  (S) 

k=l  t=0 


The  pole  locations  are  then  reflected  through  the  unit  circle  resulting  in  the  signal 
poles  within,  noise  poles  outside.  This  final  presentation  of  the  poles  describes  the 
actual  scatterer:  a  stable,  causal  system. 

For  M  segments  of  input-output  data  related  by  Equation  (8),  where 
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M  >  N  +  L  -  I,  the  relationship  is  expanded  to  matrix  form 


where  both  parts  of  the  y  \x  data  matrix  are  Toeplitz  matrices  in  which  each  row 
introduces  one  new  value  and  shifts  all  others  one  place  to  the  right.  Rewriting  this 
equation  in  more  compact  notation  yields 

D,.C^  =  3’.  W 

Extraction  is  accomplished  by  solving  this  matrix-vector  equation  for  the  ARMA 


parameter  vector, 


where  is  the  minimum  Euclidean  norm  solution  vector  for  and  is  the 
psuedoinverse  of  the  data  matrix  The  solution  vector,  Q*'*',  is  the  minimum 
length  least-squares  solution  to  the  matrix-vector  equation.  [Refs.  4  and  9] 

3.  Singular  Value  Decomposition 

The  above  solution  relies  on  singular  value  decomposition  (SVD)  of  the 
data  matrix,  to  derive  its  psuedoinverse,  This  process  separates  into  the 
product  of  three  matrices, 


where  U  and  P"are  both  orthogonal  matrices  (that  is,  lFU=r)  and  where  denotes 
the  transpose  of  V.  The  columns  of  the  matrices  1/  (of  size  Af  by  M)  and  F  (AT  by  AO 
are  composed  of  the  eigenvectors  of  D^PyJ  and  DyjDy^  respectively.  The  singular 
values  of  the  data  matrix  are  the  nonzero  components  of  the  M  by  N  diagonal  matrix, 
S.  These  singular  values  are  the  nonzero  eigenvalues  of  both  DyjDyJ  and  DyjDy^. 
Rearranging  these  components  of  the  SVD  of  Dy^  produces  its  psuedoinverse, 

2>’*'  =  (13) 

yx 

Here,  the  diagonal  components  of  S’*"  are  the  inverse  of  the  original  singular  values 
[Ref.  9]. 
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The  original  singular  values  must  be  compensatedTor  noise  if  it  is  present 
in  the  input/output  data  which  comprises  the  data  matrix.  This  compensation 
assumes  the  input  and  output  noise  variances  to  be  equal.  More  detail  on  bias 
compensation  is  provided  in  [Refs.  4  and  9]. 

B.  PERFORMANCE 

The  Cadzow-Solomon  algorithm  has  demonstrated  the  ability  to  extract  poles 
from  scattering  data.  Larison  [Ref.  9]  applied  the  algorithm  to  synthetic  target  data, 
both  measured  and  analytically  predicted  thin-wire  scattering  data,  and  aircraft  model 
data.  Murphy  [Ref.  6]  refined  the  process  by  further-application  to  synthetic  data 
sets  as  well  as  measured  thin-wire  and  aircraft  data  as  presented  in  Chapter  II  of  this 
thesis. 

1.  Noise  Tolerance 

Murphy  tested  the  synthetic  scattering  at  various  signals-to-noise  ratios 
down  to  seven  dB.  Figure  17  shows  that  the  algorithm  remained  relatively  robust,  as 
determined  by  observing  the  tightness  of  the  clustering  of  the  extracted  poles,  down 
to  20  dB  [Ref.  6].  Below  this  level  the  poles  would  not  cluster  significantly  to  lend 
a  high  confidence  to  their  location.  Figure  17  presents  fourteen  extraction  trials  of 
medium  Q  data,  each  containing  a  different  20  dB  noise  sequence  .  Notice  that  the 
higher  order  poles  (near  ttU)  are  not  extracted.  Fewer  poles  are  available  with 
successively  greater  noise  levels  [Ref.  6]. 


2.  Accuracy  Versus  Precision 

The  adequacy  of  the  extracted  poles  can  be  couched! in  terms  of  accuracy, 
the  distance  of  the  extracted  poles  to  their  true  location;  and  precisipn,  the  tightness 
of  a  group  of  poles,  regardless  of  its  location.  Figure  17  demonstrates  some 
degradation  in  accuracy  as  the  extracted  poles  are  slightly  biased  toward  the  unit 
circle  [Ref.  6]. 

Themoiseless  TDIE  data  yielded  excellent  results  in  producing  very  tightly 
grouped  poles.  This  was  to  be  expected  because  of  the  noise-free  environment  and 
the  high  quality  resonator  represented  by  the  modeled  thin  wire. 

For  actual  scattering  data,  accuracy  can  only  be  determined  through 
evaluating  the  targetdiscrimination  effectiveness  of  the  annihilation  filters  constructed 
from  the  extracted  poles.  Figure  18  shows  a  scatter  plot  of  extracted  poles  from  the 
measured  thin-wire  sca.tterer.  The  poles  are  represented  by  a  scatter  because  their 
position  can  be  perturbed  by  noise  and  is  highly  sensitive  to  the  parameters  used  in 
the  extraction  algorithm.  As  such,  the  extracted  poles  from  actual  targets  will 
represent  some  distribution  about  their  true  (unknown)  location.  Murphy  used  an 
averaging  scheme  on  the  extracted  poles  which  lie  in  tight  groups,  after  discarding 
outlying  positions,  for  a  given  angle  of  incidence.  These  averages  were  then  averaged 
over  the  available  angles  of  incidence  in  an  attempt  to  ameliorate  any  possible  bias 
due  to  the  parameters  and  to  compensate  for  the  changing  relative  excitation  of  one 
mode  over  another  with  incidence  angle.  [Ref.  6] 
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Some  indication  of  the  task  before  the  algorithm  can  be  seen  be  examining  the 
frequency  spectrum  of  that  portion  of  the  thin-wire  datarecordsused  by  the  Cadzow- 
Solomon  data  matrix.  Figure  19  shows  these  spectra  for  all  three  aspects  available. 
Notice  the  varying  levels  of  excitation  of  the  modes.  Notice  also  that  the  spectrum 
is  rather  complicated  considering  one  expects  frequencies  based  on  just  integer 
multiples  of  the  half  wavelength  matched  to  the  wire  length:  1.5,  3.0,  4.5,  6.0  GHz, 
etc.  This  is  due,  in  part,  to  the  fact  that  some  ten  data  points  at  the  end  of  early- 
time  are  included  to  gain  the  most  excited  resonances.  This  early-time  data  is  used 
by  the  algorithm  to  model  the  response  of  the  feed-forward  part  of  the  target  transfer 
function  to  the  excitation  waveform. 

3.  Aircraft  Extracted  Poles 

The  poles  extracted  from  the  model  aircraft  scattering  data  using  the 
Cadzow-Solomon  algorithm  were  much  more  highly  sensitive  to  the  input  parameters 
and  were  less  likely  to  group  as  well  as  the  thin-wire  data.  The  pole  locations,  used 
later  in  filtering,  represent  an  average  of  the  exiracted  location  for  various 
combinations  of  algorithm  parameters  and  incidence  angles.  Individual  extractions 
resulting  in  an  outlying  location  from  a  cluster  or  outside  the  unit  circle  were  ignored. 
The  highly  complex  nature  of  the  aircraft  scatterer  made  difOcult  work  for  the 
algorithm.  Figure  20  shows  the  frequency  response  of  the  portion  of  the  data  record 
processed  for  aircraft  target  two  for  the  nose-on,  tail-on,  and  thirty  degree  aspect. 
The  spectrum  is  rich  with  closely  spaced,  strong  frequency  components  at  low 
frequencies  which  the  Cadzow-Solomon  algorithm  must  sort  out.  Figure  21  shows 
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Figure  18.  Extracted  poles,  measured  thin-wire  data  [Ref.  6] 
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Figure  19.  Spectra  of  measured  thin-wire  data  processed  for  pole  extraction 
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this  same  spectral  information  for  target  three.  Notice  that  the  low  frequency  content 
of  both  aircraft  are  similarly  populated  and  complicated.  This  made  the 
extraction/filtering  discrimination  problem  very  difficult,  as  all  aircraft  models  tended 
to  have  a  similar  frequency  spectrum,  in  terms  of  many  significant  components  at 
lower  frequencies.  These  lowest  order  poles  would  ideally  be  the  ones  used  for 
target  discrimination  later  in  the  filtering  process  because  they  have  the  most  energy. 
[Ref.  6] 

C.  IMPLEMENTATION  CONSIDERATIONS 

Actual  implementation  of  a  transient  pole  extraction  facility  must  address  the 
questions  of  frequency-subject  matching  and  subject  accuracy,  signal-to-noise  ratio 
required,  and  post-extraction  pole  position  processing. 

1.  Scattering  Subject  and  Frequency 

The  scattering  objects  used  for  pole  extraction  must  faithfully  represent  the 
actual  targets  of  interest  that  will  be  encountered  by  the  radar/annihilation  filter 
system.  This  may  only  require  more  detailed  scaled  models  than  used  above  or  may 
require  full-scale  targets.  The  two  advantages  of  using  scale  models  are  the  linear 
increase  in  frequency  used  (and  linear  decrease  in  size  of  antennas  and  components 
required)  and  the  "reduced"  availability  of  full  scale  targets  —  especially  of  aircraft  of 
a  potential  adversaiy!  However,  scale  models  used  must  accurately  represent  the 
current  distribution  properties  of  the  full  scale  targets.  Validation  of  the  modeling 
may  be  conducted  on  an  aircraft  available  for  both  full-scale  and  model  testing. 
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Figure  20.  Spectra  of  target  No.  2  data  processed  for  pole  extraction 


Referring  to  the  dimensions  of  the  aircraft  targets  used  in  this.thesis,  Table 
4,  Chapter  II,  indicates  the  low  order  resonant  frequencies  of  fulf  scale  targets  lie  in 
the  7  to  70  MHz  region.  This  frequency  range  would  be  much  more  difficult  to  work 
in,  not  only  due  to  the  size  of  the  antenna  components,  but  also  due  to  the 
requirement  to  create  a  broadband  excitation  pulse  or  swept  frequency  CW  system 
in  order  to  excite  a  number  of  the  target’s  resonances.  The  ability  to  create  a  high 
power,  short  pulse  to  initiate  the  resonances  is  difficult  at  these  low  frequencies.  The 
use  of  steady-state,  stepped-frequency  excitation  may  be  a  way  to  avoid  this  problem 
[Ref.  25].  Here,  the  target  is  excited  by  relatively  long  pulses  or  CW  with 
narrow  bandwidth  subpulses,  each  centered  on  various  frequencies  in  the  resonance 
region.  The  scattered  response  is  received  coherently  for  each  frequency  and 
extruded  from  the  fbced  clutter,  transmitter,  and  receiver  response  by  deconvolution, 
in  a  similar  manner  to  that  done  in  the  TESL  (Chapter  II,  Section  B).  The  target’s 
broadband  transient  response  is  synthesized  from  the  steady-state  response  at  the 
various  transmit  frequencies  by  smoothing-in  the  frequency  domain  and  subsequently 
transforming  back  to  the  time  domain.  This  also  allows  tailoring  the  bandwidth  of 
the  excitation  to  those  frequencies  of  most  interest,  vice  transmitting  an  indiscriminate 
broadband  pulse  [Ref.  26]. 

2.  Signal-to-Noise  Ratio 

As  stated  above,  the  pole  extraction  process  requires  approximately  20  dB 
signal-to-noise  ratio.  The  extraction  process  inherently  requires  higher  SNR  than 
more  conventional  RF  imaging  from  ramp  responses  because  of  reduced  integration 
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opportunity;  the  necessity  to  process  mostly  late-time,  undriven  responses;  and  the 
sensitivity  of  the  system  poles  to  perturbations  [Ref.  27].  This  requirement 
may  dictate  the  choice  between  full  or  scale  model  targets.  The  frequency  range  of 
resonances  of  full  scale  targets  of  interest,  generally  2  MHz  to  1  GHz  [Ref.  25], 
contains  considerable  man-made  and  atmospheric  noise  [Ref.  28:pp.  34.2- 
34.13],  driving  received  noise  levels  well  above  the  predominantly  thermal  noise 
region  encountered  for  the  resonant  frequencies  of  scale  model  targets. 

3.  Post-Extraction  Position  Processing 

As  demonstrated  with  the  model  aiircraft  data,  the  extracted  pole  positions 
for  a  given  set  of  data  vary  with  the  processing  parameters  chosen.  Additionally, 
although  the  transfer  function  of  the  target  is  aspect  independent,  the  relative 
excitation  of  the  target’s  modes  affects  the  extraction  algorithm’s  output.  The  poles 
extracted  by  Murphy  [Ref.  6],  as  discussed-  in  the  previous  section,  were  gathered 
through  a  process  of  empirically  optimizing  a  set  of  extraction  parameters,  then 
dithering  them  and  averaging.  The  extracted  pole  positions  are  defined  by  some 
probabilistic  distribution,  therefore  the  "best"  set  of  poles  for  a  given  target,  relative 
to  its  actual,  unknown  poles,  is  not  defined.  Artificial  intelligence  or  an  expert  system 
would  be  a  logical  choice  to  process  and  refine  the  position  of  extracted  poles  for  a 
real  world  implementation.  The  expert  system  stage  strapped  to  the  output  of  the 
pole  extraction  algorithm  could  perform  the  sifting  of  extracted  poles  for  outliers 
(performed  by  human  eye  for  the  poles  above  [Ref.  6])  and  subsequent  averaging. 
Possibly,  a  feedback  path  could  be  initiated  to  allow  the  smart  processor  to  guide  the 
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selection  of  extraction  algorithm  parameters  to  optimize  the  extracted  positions. 
Additionally,  this  path  could  include  the  evaluation  of  the  performance  of  the 
resultant  annihilation  filters  to  refine  the  poles. 
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IV.  ANNIHILATION  FILTERING 

After  extracting  the  necessaiy  poles  for  the  targets  of  interest,  using  the 
algorithm  of  the  previous  chapter,  these  poles  were  incorporated  into  the  design  of 
annihilation  filters.  The  two  alternatives  for  this  process  explored  here  are  the  K- 
Pulse  resonance  annihilation  filter  and  the  inverse  ARMA  filter. 

A.  RESONANCE  ANNIHILATION  FILTERING 

The  scattered  response  of  a  target,  as  discussed  in  Chapter  I,  consists  of  the 
early  and  late-time  portions  as  well  as  any  noise  collected.  Morgan  describes  the 
transition  point,  T^,  with  unit  step  functions  [Ref.  7], 

y(.t)  =  (14) 

The  early-time  component  is  composed  of  excitation-dependent  physical  optics 
scattering  and  natural  mode  startups,  while  the  late-time  portion  consists  of  the 
natural  mode  resonances.  This  late-time  portion  can  be  written  as  a  sum  of  damped 
sinusoids, 

expCa^O  cos(co^r+(pp  (15) 

p=i 

where  the  amplitudes  and  phases  (complex  residues)  of  each  mode,Apexp(j(pp),  are 
aspect-dependent.  The  additive  noise  and  clutter  component,  N(t),  is  composed  of 
atmospheric  and  thermal  noise  constituents,  as  well  as  other  possible  signal  pollutants. 
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The  resonance  annihilation  filter  (RAF)  concept  employs  a  bank  of  filters, 
Figure  22,  with  the  output  of  each  being  the  result  of  convolution  of  the  filter  impulse 
response,  /i,„,  and  the  scattered  response,  y(/),  from  Equation  (14)  above. 

The  decision  criterion  for  deciding  the  type  of  target  present  is  to  match  the  filter 
output  with  the  lowest  average  late-time  energy, 

/«’»«)*  (17) 

e  =  ^ - 

“  TV-Tr 

with  the  resonances  embodied  in  that  filter.  Here  is  the  start  of  late-time,  picked 
late  enough  to  ensure  no  early-time  interaction  remains,  and  is  the  end  of  the 
integration  time,  which  would  ideally  approach  infinity  for  maximum  discrimination, 
but  must  be  time-limited  to  prevent  noise  from  obliterating  the  resonant  energy.  This 
late-time  start  point  is  extended  beyond  the  target  late-time  discussed  in  Chapter  I, 
Section  A,  by  the  delay  of  the  annihilation  filter,  Here  the  input  signal,  y(t),  is 
normalized  by  its  late-time  RMS  energy  and  the  filter  response,  /i,„,  is  normalized  by 
its  RMS  value  over  its  record  length  [Ref.  10].  Morgan  discusses  the  requisite 
features  of  the  filters,  li,„,  in  [Ref.7]  and  points  out  that  using  a  finite  impulse 
response  (FIR)  filter  precludes  the  possibility  of  early-time  energy  being  carried  over 
into  late-time  and  polluting  the  target  discrimination  process. 


51 


B.  THE  K-PULSE  FILTER 


1.  Methodology 

In  order  to  realize  the  annihilation  filter  processing  with  digital  computers, 
the  discrete  time  filter  impulse  response,  must  be  solved  for.  The  following 
derivation  is  due  to  Morgan  [Ref.  29].  First  the  late-time  response  in 
Equation  (15)  is  written  in  discrete  time  format,  utilizing  the  complex  z-transform 
operator  to  represent  the  sinusoidal  poles, 

p^N 

po-N 

p*0 

where  2p=exp(5pAr),  z_p=z*  =exp(5*  A/),  s  represents  each  pole  in  the  Laplace 
domain,  and  is  the  combination  of  the  magnitude  and  residue  terms  which  scale 
each  pole’s  relative  contribution  to  the  scattering.  Additional  flexibility  can  be 
incorporated  by  allowing  k  integer  delays  between  filter  weights,  after  which  the 
convolution  process  of  Equation  (16)  is  written  as  a  difference  equation 

w(/i)  =  £  (1^) 

After  substituting  the  late-time  response  of  Equation  (18)  and  rearranging  the  order 
of  summation,  the  equation  below  is  set  to  zero,  indicating  the  annihilation  action  on 
the  late-time  energy. 
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P’N  /m»M  \ 

M")  =  E  Vp"  E 

p=-jv  ) 

pi‘0 


=  0. 


(20) 


The  coefficient,  Iiq,  can  be  set  to  normalize  the  noise  power  passed  through  the  filter, 
yielding  the  result 


m=M 


E  =  K- 


(21) 


mcl 


This  result  can  be  adjusted  to  represent  a  real  linear  system  while  reducing  the 
dynamic  range  of  the  coefficients,  which  results  in  two  sets  of  N  equations 


E {exp(a^Af)p2 


m=M 

E 

m=l 


M 

U 


—TO 


I. 

=  -hn  cosw  —  0) 

0  n  2 ' 


Ar  (22a) 


E  K  {exp(o^Af)}*(2  )  sin|s|y-TOjw^Af|  =  -h^  sin^  j 


where  p  =  1  io  N.  To  form  a  square  matrix  of  equations  which  results  in  exact 
cancellation,  M  must  equal  2N.  This  derivation  essentially  places  filter  zeros  upon 
the  poles  of  the  input  signal  to  be  annihilated. 

Simon  [Ref.  10]  discusses  a  different  approach  to  achieve  an  equivalent 
result.  The  coefficients,  Op  and  /3p,  are  solved  for  each  mode  by  substituting  the 
discrete  time  version  of  Equation  (14)  into  a  three-point  difference  equation 
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fl^pyi[(«-l)Af]+yJnAr]+/?p»J(«+l)Af]  =  0.  (23) 

This  result  yields  an  individual  three  point  FIR  filter  for  each  mode.  Thez-transform 
derivation  outlined  above  is  equivalent  to  combining  these  in  cascade  and  overlapping 
the  leading  and  trailing  points  for  each  filter  to  form  a  single  filter  to  handle  all  the 
target’s  modes. 

The  program  to  create  the  K-Pulse  filter  impulse  response  based  on  pole 
pair  input  is  in  Appendix  D.  Figure  23  contains  an  example  K-Pulse  filter  impulse 
response  created  for  the  medium  Q  synthetic  pole  set;  the  other  K-Pulse  filters  look 
similar. 

2.  Smoothing 

The  K-Pulse  filter  algorithm  is  difference  equation  based,  which  results  in 
noise  amplification  through  the  filter.  Notice  that  in  Figure  24,  which  displays  the 
spectrum  of  the  filter  of  Figure  23,  that  the  filter  rises  with  frequency  past  the  region 
where  it  performs  the  annihilation.  In  fact,  the  dynamic  range  of  this.plot  precludes 
examining  the  structure  of  the  filter’s  spectral  magnitude  at  low  frequencies.  The 
frequenc}'  rise  at  the  high  frequency  end  of  the  spectrum,  after  the  last  zero,  goes  as 
where  N  is  the  number  of  pole  pair  zeros  incorporated  in  the  filter  [Ref.  29].  If 
the  sampling  bandwidth  cannot  be  made  small  enough  that  the  last  zero  is  near  the 
fJ2  folding  frequency,  smoothing  of  the  filter’s  impulse  response  must  be 
accomplished  before  it  can  be  used  on  data  with  noise.  Simon  convolved  the  K-Pulse 
filter  impulse  response  with  a  double  Gaussian  smoothing  function  [Ref.  10). 
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Magnitude 


Figure  24.  Spectrum  of  medium  Q  K-Pulse  filter 
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Referring  back  to  Figure  4,  Chapter  II,  the  spectrum  of  the  double  Gaussian  is  band- 
limited  with  the  high  and  low  frequency  roll-offs  determined  by  the  width.of  the  rapid 
and  narrow  Gaussian  pulses  used  in  its  formation.  The  high  frequency  tail  drops  off 
as  a  Gaussian,  that  is  exp(-a/),  where  or  is  a  coefficient  determined  by  the  relative 
widths  of  the  pulses.  This  decaying  exponential  is  used  to  overcome  the  polynomial 
growth  of  the  K-Pulse  filter.  The  programs  used  to  create  the  smoothing  functions 
and  the  smoothed  K-Pulse  filter  are  in  Appendices  E  and  F.  The  bandwidth  of  the 
smoothing  functions  used  was  determined  by  qualitatively  examining  the  spectrum  of 
the  smoothed  pulse.  The  desired  filter  frequency  response  would  indicate  the  highest 
frequency  zero  with  a  definitive  local  minimum  and  then  the  spectrum  would  rapidly 
fall  off,  driven  by  the  double  Gaussian.  The  smoothed  version  of  the  ARMA  medium 
Q  synthetic  data  K-Pulse  of  Figure  23  is  in  Figure  25;  its  spectrum  is  Figure  26. 
Smoothed  K-Pulse  filters  were  used  to  test  the  annihilation  filter  concept  on  all  the 
various  signals  discussed  in  Chapter  II. 

3.  Performance 

a.  ARMA  Synthetic  Scattering 

To  evaluate  the  effectiveness  of  the  K-Pulse  filter,  scattering  data  was 
created  for  the  poles  sets  indicated  in  Chapter  II,  as  well  as  for  poles  both  five 
percent  above  and  below  the  baseline  sets,  at  infinite  (no  noise),  20, 15, 10  and  7  dB 
signal-to-noise  ratio.  These  data  records  were  convolved  with  the  K-Pulse  filters 
created  for  the  "correct"  pole  pairs  using  the  convolution  routine  in  Appendix  G.  As 
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Figure  26.  Spectrum  of  smoothed  medium  Q  K-Pulse  filter 
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a  representative  case,  Figure  27  shows  the  output  of  the  medium  Q  smoothed  K- 
Pulse  filter  operating  on  the  20  dB  correct  pole  scattering  data.  Figure  28:shows  the 
same  filter  processing  the  20  dB  five  percent  above  data  and  Figure  29  shows  the 
same  filter  processing  the  correct  data  at  7  dB  SNR.  Table  6  shows  the  resulting 
ratios  of  average  late-time  energies  between  the  correct  and  incorrect  target/filter 
matching. 

TABLE  6.  RATIOS  OF  LATE-TIME  ENERGIES:  ARMA  SYNTHETIC  DATA 


Test  Signal 

Signal  to 
Filter 
Matcli 

SNR 

00 

20  dB 

15  dB 

10  dB 

7dB 

High  Q' 

5  %  above 

28243 

11.03 

4.25 

2.07 

1.56 

5  %  below 

74442 

26.80 

8.53 

3.59 

2.29 

Medium 

5  %  above. 

20639 

4.05 

1.91 

1.26 

1.03- 

5  %  below 

67721 

11.79 

4.4 

2.07 

1.53 

Lx)wQ*’ 

5  %  above 

3568 

0,98 

0,98 

0.99 

0.99 

5  %  below 

8814 

1.03 

0.99 

0.99 

0.99 

Note  1:  Integration  time  =  5.14  ns 
Note  2:  Integration  time  =  1.64  ns 


Late-time  commences  at  r^=2.86  nanoseconds  for  these  trials.  This 
includes  the  0.3  nanosecond  incident  pulse  width  after  a  0.2  nanosecond  delay, 
approximately  0.78  nanoseconds  resulting  from  20  stages  in  the  ARMA  signal 
generator  which  represents  the  2D/c  propagation  delay  over  the  scatterer,  and  an 
additional  1.57  nanoseconds  of  "shadow-time"  created  by  the  delay  of  the  smoothed 
K-Pulse  filter.  The  length  of  the  integration  times  beyond  the  start  of  late-time  was 
determined  by  observing  the  result  of  the  mismatched  filter.  The  summation  period 
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was  terminated  at  a  point  just  past  where  the  significant  signal  structure  remained  in 
the  output  of  an  incorrectly  matched  filter  to  avoid  noise  "swamping  out" 
discrimination.  This  was  determined  to  be  eight  nanoseconds  into  the  time  record 
for  the  high  Q  data  and  4.5  nanoseconds  for  the  medium  and  low  Q  case;  this 
resulted  in  integration  periods  of  5.d4  and  1.64  nanoseconds,  respectively.  Ratios  of 
3  dB  are  still  achievable  down  to  10  dB  SNR  for  the  medium  Q  five  percent  below 
target.  This  indicates  the  validity  ofithe  K-Pulse  filtering  discrimination  for  a  simple 
target  case  in  the  presence  of  noise. 

b.  TDIE  Thin-Wire  Scattering 

The  synthetic  thin-wire  data  was  evaluated  by  creating  a  K-Pulse  filter 
matched  to  the  extracted  poles,  as  well  as  filters  derived  from  five  percent  variations 
in  the  extracted  poles.  Because  of  the  lower  inherent  sampling  frequency  of  the 
TDIE  program  which  created  these  signals,  480  samples  in  20  nanoseconds,  and  the 
fact  that  poles  were  extracted  at  sufficiently  high  enough  frequencies,  these  kill  pulses 
did  not  require  smoothing.  The  correct  filters  were  convolved  with  both  correct  and 
incorrect  data  at  various  signal-to-noise  ratios.  The  results  are  summarized  in  Table 
7,  The  late-time  integration  ran  between  1.0  and  3.0  nanoseconds. 

The  results  of  the  thin-wire  data  point  to  the  aspect  independence  of 
the  annihilation  filter  concept.  Notice  that  the  ratios  are  greatest  for  low  aspect 
angles.  This  is  due  to  the  fact  that  a  single  late-time  start  value  was  used,  when  in 
fact,  for  the  broadside  (90  degree  case)  there  is  essentially  no  2D/c  contribution  to 
late-time.  Therefore,  the  latest  late-time  value  selected,  which  is  based  on  the  longest 
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Figure  27.  Filter  output  for  correct,  medium  Q  synthetic  data,  20  dB  SNR 
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Figure  28.  Filter  output  for  five  percent  above,  medium  Q  data, 
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TABLE  7.  RATIOS  OF  LATE-TIME  ENERGIES:  TDIE  THIN-WIRE  DATA 


Target  Aspect 
(Degrees) 

Filter  to 
Target 
Match 

SNR 

00 

20 

15 

10 

7 

5 

S%  above 

214.58 

11.72 

1.59 

Z16 

1.58 

5%  below 

462.50 

24.46 

3.10 

2.42 

30 

5%  above 

50.92 

1Z52 

5.68 

Z73 

1.94 

5%  below 

96.37 

21.09 

8.68 

3.60 

2.32 

45 

5%  above 

91.80 

14.92 

6.44 

2.92 

2.00 

S%  below 

132.90 

19.94 

8.04 

3.30 

2.13 

70 

5%  above 

5.31 

4.56 

3.60 

2.39 

1.82 

5%  below 

5.33 

4.72 

3.82 

2.64 

2.06 

90 

5%  above 

56.59 

3.36 

1.82 

1.27 

1.13 

5%  below 

49.94 

2.48 

1.42 

1.12 

1.07 

target  presentation  aspect,  starts  much  later  than  the  significant  late-time  portion  for 
the  broadside  case.  Additionally,  for  the  symmetric  broadside  aspect,  only  the  odd 
modes  are  excited,  further  reducing  the  discrimination. 

Figure  30  contains  the  correct  and  five  percent  above  filter  outputs 
for  the  30  degree  aspect  case  at  15  dB  SNR.  Notice  the  residual  energy  past  one 
nanosecond  for  the  incorrectly  matched  filter.  Figure  31  shows  the  spectra  for  the 
late-time  of  these  two  output  data  records. 

Measured  Thin-Wire  Scattering 

The  same  five  percent  variation  of  all-aspect,  averaged  poles  was 
incorporated  into  filters  created  for  the  measured  thin-wire  data.  The  input  data 
records  for  all  three  aspects  (30,  45  and  90  degrees)  were  passed  to  all  three  filters. 


66 


\v(t)  (Volts) 


Figure  30.  TDIE  thin-wire  filter  output,  30° 
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The  resulting  ratios  of  late-time  energies  are  in  Table  8. 

The  best  case  for  this  data  was  30  degrees  aspect,  comparing  the  five 
percent  below  to  the  correct  filter.  The  comparisons  of  the  annihilation  filter  outputs 
in  the  time  and  frequency  domains  are  in  Figures  32  and  33,  respectively.  Notice  that 
even  the  correct  filter  output  still  contains  some  frequency  information.  This  is  due 
to  the  biasing  of  the  extracted  poles  from  their  true  location  in  the  Cadzow-Solomon 
extraction  process.  These  relatively  small  inaccuracies,  however,  do  not  preclude 
proper  identification  of  members  of  this  simple  class  of  targets. 

TABLE  8.  RATIOS  OF  LATE-TIME  ENERGIES:  MEASURED  THIN-WIRE 


Tai;get  Asptct  (Degrees) 

Filter  to  Target  Match 

Ratio 

30 

5%  above 

10.22 

5%  bciow 

1Z6S 

45 

5%  above 

4.81 

5%  below 

636 

90 

5%  above 

Z71 

5%  below 

Z23 

d.  Aircraft  Scattering 

The  extracted  poles  for  aircraft  target  No.  2  were  incorporated  into 
an  annihilation  filter,  which  was  then  smoothed.  Scattering  signals  for  various  aspect 
angles  of  targets  No.  2  and  3  were  passed  through  the  filter.  Late-time  energj' 
summation  was  performed  for  ^vo  nanoseconds  starting  at  the  appropriate  time, 
based  on  the  late-time  inherent  in  the  scattering  signal  plus  a  processing  "shadow¬ 
time"  resulting  from  the  K-Pulse  filter.  The  smoothed  K-Pulse  filter  for  target  No. 
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2  is  shown  in  Figure  34,  its  spectrum,  in  Figure  35.  The  results  of  this  test,  in  Table 
9,  show  that  the  filter  could  not  correctly  discriminate  target  No.  2  from  target  No. 
3.  This  problem  again  relates  to  the  difficulty  in  extracting  all  the  low  frequency 
poles  for  each  target  and  the  similar  high  population  of  significant  frequency 
components  in  this  region  for  both  targets,  as  discussed  in  Chapter  III.  Figure  36 
shows  the  output  waveforms  for  targets  No.  2  and  No.  3,  each  with  nose-on  incidence 
over  the  late-time  energy  integration  period.  Figure  37  shows  the  spectra  of  these 
output  data  records.  Notice  that  a  pole  at  approximately  2.5  GHz  was  not 
annihilated  in  the  correct  target  signal.  This  was  due  to  the  fact  that  this  pole  was 
extracted  for  both  targets,  and,  as  such,  was  not  included  in  the  target  No.  2  filter. 
At  this  specific  aspect,  however,  it  is  excited  in  target  No.  2  and  not  target  No.  3, 
which  leads  to  incorrect  identification.  Figures  38  and  39  contain  the  time  and 
frequency  domain  plots  for  the  best  case  indie  ited  in  Table  9,  the  planform  view  with 
the  wings  parallel  to  the  electric  field.  Here,  while  the  2.5  GHz  resonance  seems 
present  in  both  aircraft,  it  does  not  preclude  correct  identification.  These  late-time 
only  FFT  plots  indicate  that  the  effective  signal-to-noise  ratio  in  late-time  is  much  less 
than  that  if  the  entire  signal,  including  the  strong  early  time  return,  is  included. 
Generally,  for  all  the  targets  examined  in  this  thesis,  the  ratio  of  early  to  late-time 
peak  power  is  approximately  12  to  16  dB,  therefore  the  late-time  signal-to-noise  ratio 
is  less  than  the  early-time  SNR  by  that  much. 

The  observations  above  relate  to  another  result  encountered  here, 
namely,  a  lack  of  frequency  re.solution  in  the  K-Pulse  filter  algorithm  was  detected. 
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Frequency  (GHz) 


Figure  35.  Spectrum  of  smoothed  K-Pulse  filter,  target  No.  2 


74 


Time  (nanoseconds) 


Figure  38.  Late-time  output  of  targets  No.  2  and  No.  3  from  target  No.  2 
filter,  wings  parallel 


Figure  39,  Spectra  of  output  of  targets  No.  2  and  No.  3  from  target  No.  2 
filter,  wings  parallel 
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TABLE  9.  RATIOS  OF  LATE-TIME  ENERGIES:  AIRCRAFT  DATA 


Target  Aspect  (Degrees) 

Ratio  (Tgt  3  LTE/Tgt  2  LTE) 

Nose  On 

0.20 

30 

0.38 

90 

0.21 

Tail  On 

3.33 

Pianform,  Wings  Parallel 

3.38 

Planform,  Fuselage  Parallel 

1.32 

These  targets  (targets  No.  2  and  No.  3)  present  very  closely  spaced  poles  upon 
extraction,  as  discussed  in  Chapter  III,  Section  B.  The  target  No.  2  filter  used  for  the 
filtering  results  in  Table  8  contains  a  zero  at  1.3  and  1.5  GHz.  Examination  of  a 
"zoomed- in"  portion  of  the  FFT  of  the  unsmoothed  version  of  this  filter.  Figure  40, 
shows  that  the  zero  corresponding  to  1,5  GHz  is  not  clearly  distinguishable  from  that 
at  1.3  GHz.  This  means  that  any  target  with  a  pole  in  this  neighborhood  will  have 
its  energy  annihilated,  thus  reducing  the  probability  of  correct  discrimination.  For 
this  reason,  the  pole  at  approximately  2.5  GHz,  separated  between  targets  two  and 
three  by  less  than  the  poles  discussed  here,  precluded  incorporation  of  this  zero  into 
the  target  two  filter. 

e.  Summary 

These  results  are  not  encouraging,  however,  further  work  in  the 
refinement  of  the  extraction  process  and  improvements  in  the  annihilation  filters  may 
enable  the  use  of  this  concept  on  real  scattering  data,  either  from  models  or  full-scale 
aircraft,  down  to  near  the  signal-to-noise  ratios  realized  by  the  synthetic  data.  The 
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Figure  40.  Close-up  of  FFT  of  target  No.  2  K-Pulse  filter 
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aircraft  targets’  frequency  responses  reveal  them  to  be  extremely  complicated 
scatterers  relative  to  the  other,  simpler  models  used  in  this  thesis.  Concentration  on 
intermediate  levels  of  target  complexity  will  allow  these  algorithms  to  mature  or 
other,  more  applicable,  ones  to  surface.  One  possible  avenue  of  investigation  is 
reducing  the  bandwidth  of  the  excitation  in  order  to  limit  the  number  of  similar  poles 
extracted  and  filtered  for  each  different  target.  Investigation  of  the  frequency  content 
of  scatterers  in  smaller  bandwidths  should  be  conducted  to  discern  if  there  is  an 
optimum  frequency  range  where  various  target  frequency  responses  offer  the 
maximum  potential  for  discrimination. 

4.  Implementation  Considerations 

The  relevant  concerns  for  application  of  the  resonance  annihilation  filter 
radar  target  classification  method  in  an  aptual  system  include  frequency,  bandwidth, 
SNR  required  and  its  affect  on  transmit  power,  antenna  design,  clutter, 
electromagnetic  interference,  and  post-filter  processing  and  display. 

a.  Frequency 

The  transmit  frequency  for  the  annihilation  filter  radar  must  coincide 
with  the  resonances  of  full-scale  targets  -  approximately  7  to  70  MHz  for  fighter-sized 
aircraft.  The  same  considerations  raised  earlier  for  the  scattering  facility  apply  here 
(Chapter  III,  Section  C),  except  that  the  option  to  scale  the  subjects  to  microwave 
frequencies  no  longer  exists.  The  difficulties  encountered  at  this  frequency  range 
include  increased  noise  levels,  the  requirement  for  large  antennas  to  create  high  gain, 


a  crowded  spectrum  of  commuriications  users,  and  difficult  propagation  properties 
as  the  bandwidth  spans  the  middle  HF  region,  which  includes  surface  and  sky-wave 
modes  of  propagation,  to  the  lower  VHF  region,  which  is  liiie-of-sight 
[Ref.  30,  p.  529]. 

b>  Bandwidth 

In  order  to  effect  discrimination  of  similarly  sized  targets,  presumably 
with  only  slightly  different  pole  locations  within  a  band,  the  system  bandwidth  must 
be  sufficient  to  excite  some  number  of  poles,  possibly  two  to  four.  As  stated  above, 
the  frequencies  of  interest  cover  a  60  MHz  bandwidth.  Auton  indicates  an  even 
wider  bandwidth  of  20  MHz  to  1  GHz  and  suggests  a  stepped  frequency  system  which 
is  the  most  economical  way  of  creating  these  wide  bandwidths  in  this  frequency  range 
[Ref.  25].  However,  as  suggested  above,  perhaps  a  subset  of  this  wide  band  would 
be  an  area  in  the  frequency  domain  where  maximum  discrimination  is  possible  among 
similarly  sized  targets.  If  so,  this  would  ease  the  bandwidth  requirements  discussed 
here. 

Pole  extraction  by  Murphy  [Ref.  6]  tends  to  indicate  that  the  poles  of 
different  scatterers,  while  different  from  one  another,  tend  to  reside  in  the  same 
neighborhood  for  each  successive  mode.  While  this  tends  to  make  target 
discrimination  more  difficult,  it  could  allow  skipping  all  but  the  most  relevant 
frequencies  in  a  stepped  frequency  system  such  as  that  proposed  by  Auton  [Ref.  25]. 
This  skipping  would  save  processing  time  and  reduce  interference.  The  transformed 
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transient  time  response  could  be  "filled-in"  with  the  missing  frequent  through  fitting 
algorithms  during  the  deconvolution  process  [Ref.  25]. 

Auton’s  resonance  based  target  identification  system  specifically 
proposes  a  bank  of  parallel  CW  transmit/receive  systems  each  operating  at  successive 
frequencies  in  the  band  of  interest,  each  with  a  narrow  receiver  bandwidth  of 
approximately  0.1  Hz.  The  system  relies  on  a  companion  microwave  radar  to 
generate  range-rate  information  to  center  the  narrow  receiver  on  the  target  doppler 
and  avoid  clutter.  The  return  signal  is  processed  through  deconvolution,  as  described 
earlier  (Chapter  III,  Section  C).  Auton  then  proposes  a  rational  fitting  algorithm  to 
determine  the  pole  locations  for  comparison  to  known  locations.  Here,  the 
substitution  of  a  Resonance  Annihilation  Filter  bank  may  be  made  to  speed  the  final 
step  of  the  classification  process.  [Ref.  25] 

c.  SignaUto-Noise  Ratio  and  Range 

As  a  minimum,  a  late-time  signal-to-noise  ratio  of  15  dB  vtill  be 
assumed  for  target  identification  via  the  RAF  concept.  Auton  [Ref.  25]  cites  26  dB 
at  the  end  of  an  approximately  ten  second  integration  time.  This  has  implications  for 
the  range  equation  applied  to  this  radar.  Assuming  that  the  resonance  based 
identification  radar  will  be  used  with  an  accompanying  microwave  surveillance  radar, 
it  can  be  confined  to  operate  as  a  surface-wave  radar  in  the  HF  portion  of  the  band 
and  a  line-of-sight  VHF  radar  in  that  band.  The  range  equation  for  this  case 
(ignoring  ionospheric  losses)  is  [Ref.  30,  pp.  530-531] 
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where  R  is  range,  is  average  power,  G,  anci  are  the  f  "nsmit  and  receive 
antenna  gains.  A,  is  the  wavelength,  cr;is  target  cross  section,  7^  is  coherent  processing 
time,  k  is  Boltzmann’s  constant,  is  received  noise  level  above  kTB  due  to  other 
noise  sources  (discussed  below),  {SIN)  is  the  required  signal-to-noise  ratio,  and 
accounts  for  system  loss. 

This  frequency  range,  20  to  100  MHz,  is  dominated  ;in  the  low  end  by 
atmospheric  noise  and  at  higher  frequencies,  by  galactic  noise.  Assuming  a  rather 
large  system,  such  as  this  radar  identification  system,  woujd  belocated  in  a  rural  area 
to  avoid  man-made  noise,  the  values  of  F^,  in  decibels  above  ktB,  range  from  32  dB 
at  10  MHz  down  to  7  dB  at  90  MHz  [Ref.28:p.  34.9,].  Substituting  these  values  into 
the  range  equation  above  yields  Figure  41.  MauTmade  suburban  or  urban  noise 
levels,  if  included  in  the  above  calculation,  would  cause  the  power  levels  required  to 
be  some  20  to  30  dB  higher  for  frequencies  above  25  MHz,  The  other  values 
involved  in  the  equation  are  included  in  the  figure  and  each  is  discussed  below. 
Nominally,  such  a  system  would  require  approximately  250  kW  average  power  in 
order  to  sufficiently  the  excite  the  less  damped  poles  in  the  low  frequency  end  of  the 
resonance  spectrum  and  to  allow  for  increased  man-made  noise  levels  above  galactic 
noise  in  the  frequency  region  above  25  MHz.  More  power  may  be  required  at  the 
lower  frequencies  even  still  to  overcome  the  losses  due  to  surface-wave  attenuation 
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caused  by  the  radar  energy  coupling  into  the  ground;  This  phenomenon  is  most 
critical  around  15  MHz  and  falls  off  with  increasing  frequency  [Ref.  31]. 

The  radar  cross  section  value  used  above  represents  an  average;  in 
fact,  as  this  is  the  resonance  region  for  the  targets  of  interest,  the  return  signal 
strength  can  vary  greatly  with  small  changes  in  frequency.  As  the  radar  generates  a 
wideband  signal,  the  stronger  resonances,  used  fo&identification,  will  boost  the  return 
signal.  The  average  value  indicated  is  derived  from  averaging  a  value  at  each 
frequency  determined  by  approximating  the  actual  target  backscatter  cross  section 
with  that  of  a  2:1  prolate  spheroid  with  perpendicular  electromagnetic  excitation 
incident  at  30  degrees  [Ref.  32:p.  361]. 

d.  Antenna  Design 

The  necessity  to  operate  over  a  ten-to-one  bandwidth  will  deihand 
much  from  the  antenna  system  used  for  a  resonance  identification  radar.  The  gain 
value  used  above,  although  realistic  [Ref;  30:p.  532],  will  not  remain  constant  over  the 
band.  As  such,  the  antenna  will  represent  a  compromise  design  rather  than  a  point 
design. 

The  polarization  of  the  system  would  ideally  be  horizontal,  as  this 
would  allow  for  the  most- coupling  into  the  resonances  of  the  fuselage  and  wingspan 
of  the  aircraft  targets  in  level  flight.  However,  due  to  the  size  of  the  arrays  involved 
and  the  higher  surface-wave  attenuation  for  horizontal  polarization  for  the  lower 
frequency  range  of  the  system,  the  .most  suitable  antenna  will  likely  be  a  large  array 
of  vertical  monopoles  or  dipoles. 
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Figure  41.  Average  power  versus  frequency  and  SNR 


In  order  to.  reduce  the  clutter  and  noise,  the  radar  should  have  as 
narrow  azimuth  and  elevation  beamwidths  as  possible.  However,  the  large  antenna= 
size  will  limit  this:  a  two  degree  azimuth  beamwidth  (broadside)  wouldirequire  an; 
aperture  length  of  1700  meters  at  10  MHz  or  190  meters  at  90  MHz,  One 
alternative  often  exploited  at  these  low  frequencies  is  the  use  of  a  broad^  low  gain 
transmit  beamwidth  generated  by  driving  only  a  few  elements  of  the  array,  overlaid^ 
with  a  number  of  more  precise  receive  beamwidths.  This  reduces  the  cost  and= 
complexity  and^still  retains  some  gain  and  resolution  cell  precision  [Ref.  30:p.  531]. 

The  varying  propagation  effects  between  the  frequencies  above  and 
below  40  MHz  may  warrant  two  separate  antenna  systems-tied.together  for  wideband- 
processing  and  annihilation  filtering.  Here,  the  stepped  frequency  nature  of  the 
system  suggested  by  Auton:[Ref.  25]  may  make  this  possible.  The  higher  frequency 
resonances  of  the  target  scatterers  are  the  most  highly  damped,  so  it  may  be  found 
upon  testing  full-scale  aircraft  in  actualconditions  that  the  the:  higher  end  of  the  band 
is  not  cost  effective. 

Another  consideration  in  the  design  and  placement  of  the  antenna  is 
the  interference  phenomena  due  to  multipath  prorpagation  from  reflected  signals  at 
low  grazing  angles.  The  patterns  which  develop  are  highly  frequency  and  geometry 
dependent.  At  tlie  maximum  extent  of  the  effect,  the  return  signal  strength  can  be 
doubled  or  eliminated  [Ref.  33:pp.  341-361].  This  effect  may  be  ameliorated 
by  the  wideband  nature  of  the  system,  but  may  also  inconsistently  attenuate  critical 
resonances  required  for  target  discrimination. 


e.  Clutter 

The  problem  of  clutter  will  be  even  moresinsidiousiin  the  resonance 
identification  radar  than  in  normal  search  or  tracldng  radar  in  that  the  frequency 
content  itself  of  the  return  signal  is  ‘mportant.  Auton  suggests  the  use  of 
narrowband,  doppler  tracking  filters  to  avoid  most  of  the  clutter  [Ref.  25]. 
Additionally,  controlling  the  propagation  and  precise-timergating  may  reduce  clutter 
from  its  various  sources  over  the  band:  reflections  off  the  surface  and  sky-wave 
reflection  from  the  ionosphere  at  higher  HF  frequencies. 

Sea  clutter  is,  in  fact  lower,  in  the  frequency  region  of  interest  than 
in  the  microwave  frequencies,  however,  the  advantage  is  more  pronounced  with 
horizontal  polarization,  which  may  not  be  an  option,  considering  the  element  size. 
In  an  attempt  to  classify  low  flying  targets  Over  water,  controlling  the  sea  clutter  by 
manipulating  elevation  angle  may  become  difficult  in  light  of  the  requirement  to 
avoid  multipath  interference.  [Ref.  30:p*  483] 

The  scattering  data  collected  for  this  work  uses  a  deconvolution 
process  to  remove  non-target  returns.  This  processmay  have  application  in  the  real 
world  systems  whereby  the  system  could  alternate  the  timing  or  direction  of  the 
information  processed  to  allow  a  comparison  of  the  target  return  and  a  non-target 
return  via  deconvolution.  This  processing  would  have  to  be  done  rapidly  and  may 
complicate  the  annihilation  filtering  process,  but  may  be  a  necessary  clutter  rejection 
tactic. 
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f.  Electromagnetic  Interference 

Examination  of  a  frequency  allocation  chart  between  7  and  70:MHz 
indicates  a  crowded  spectrum  of  communications  users  filling  the  band; 
[Ref.  34:p.  1.10].  Also,  due  to  ionospheric  effects,  the  generation:  of  high; 
power,  wideband;  transmission  by  this  system  will  create  interference  problems. 
Again,  a;  stepped  frequency  system  would  allow  skipping  protected  or  overcrowded 
portions  of  the  band  with  subsequent  deconvolution  and  rational  function  fitting  filling 
in.  Additionally,  during  operations  of  this  system,  friendly  communications  may  be 
coordinated  to  avoid  interference. 

As  a  victim  of  EMI,  the  system  must  deal  with  the  crowded  spectrum 
and  man-made  noise,  as  discussed  earlier.  This  drives  up  the  required-  transmit 
power  for  a  given  SNR,  and  will  dictate  deployment  location  in  remote^  yet,  tactically 
significant  areas.  As  a  victim  of  electronic  countermeasures,  the  system’s  wide 
bandwidth  may  offer  some  inherent  jam  resistance.  However,  the  criticjility  of  the 
various  frequency  components  used  to  discriminate  one  target  from  another  and  the 
requirement  for  a  greater  SNR  than  required  for  mere  detection  systems,  may  make 
the  system  vulnerable.  Operationally,  the  loss  of  this  system  to  jamming  is,  itself,  a 
hostile  target  identifier. 

g.  Post-filter  Processing  and  Display 

The  annihilation  filters  operate  on  the  decision  criterion  of  ratios  of 
late-time  energy.  In  the  presence  of  noise,  these  ratios  are  probablistic  in  nature. 
Norton  conducted  a  rigorous  analysis  of  the  probability  density  function  of  this 
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process  in  [Ref.  4].  In  order  to  operate  as  a  real-time  tacticalipr  strategic  decisioii 
aid  the"  target  resonance  radar  must  quickly,  indicate  its  classification:  information  to 
the  operator,  yet  it  must  account  for  the  possibilities  of  incorrect  identification^  or 
inconclusive  results.  This  will  require  high  speed  processing  of  the  resonance 
annihilation  filter  outputs  followed  by  the  appropriate  display.  This  display  may  take 
the  form  of  a  ranking  of  the  most  probable  identities  of  the  target  being  tracked  and 
the  confidence  level  of  the  identifications  listed.  Pre^sent  work  in  artificial  intelligence 
systems  applied  to  imaging  systems  such  as  synthetic  aperture  radar  offer  a  lead  in 
the  area  of  real-time  man-machine  interface;  for  rapid  target  identification. 

h.  Summary 

The  implementation  of  a  natural  respnance  target  identification  system 
represents  no  greater  technological  risk  than  other  radar  systems  that  have  been 
developed  in  this  frequency  range.  The  mp_st  difficult  iniplementatidh  task  will  be 
isolating  the  relevant  target  return  for'identificatioh  processingTrom  the  entirety  of 
the  radar  return  in  time  and  space,  in  order  to  exclude  as  much  noise  and  clutter  as 
possible  and  to  concentrate  the  analysis  on.  a  single  target. 

G.  INVERSE  ARMA  FILTER 
1.  Methodology 

The  inverse  ARMA  filter  represents  an  extension  of  the  K-Pulse  filter 
concept  to  include  canceling  early-time  resonance  energy  as  well  as  late-time.  As 
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such,  it  immediately  sufrenders  thejhspect  independence  property  of  the  annihilation^ 
filter  discussed  earlier. 

The  transfer  function^ofiarscatterer  is  repre'sented  as  h  ratio  ofz  transform 
polynomials 


bn  +  b,z~^  +  •  •  •  +  6,z"^ 
1  +  ajz  *  +  •*•  + 


(25) 


I'he  inverse  ARMA  filter  frequency  response  is  formed  by  inverting  this  transfer 
function 


Hniz)  = 


+  + 


%  1 

1  +  -4.z-‘,+ 


h  .1 

•  +  ~z  ^ 


(26) 


Consider  the  scattered  response  from  the  target  y(/)=/i,*A:(/),  due  to  the  interaction 
of  the  incident  pulse,  x(t),  ^vith  the  target’s  impulse  response,  /z,.  In  the  frequency 
domain,  Y(z)=H,(z)X(z).  Now,  the  inverse  ARMA  filter  impulse  response,  /z^^,  is 
derived  from  the  inverse  of  the  target’s  transfer  function,  so  The  Altering 

action  is  described  by  zv(/)=y(z)*/z^,  Mich  gives  this  result  in  the  frequency  domain: 


IKCz)  =  Y(z)Hu(.z)  =  Ymi\z)  =  H,iz)xm;'Cz)  =  X(.z).  (2’) 

In  fact,  the  appropriately  matched  inverse  ARMA  filter  will  annihilate  all  early  and 
late-time  resonant  energy,  leaving  only  a' replica  of  the  incident  pulse.  For  an  inverse 
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ARMA  filter  which  is  not  matchedito  the  target,  ,,110  such  cancellation  will 

occur.  [Ref.  35] 

The  inverse  ARMA  filter  achieves  another,  ^advantage  over  the  K-Pulse 
filter  in  that  it  is,  itself,  an  ARMA  filter  that  inherently  smoothes  noisy  input.  As 
such,  it  requires  no  additional,  smoothing  as  the  K-Pujse  filter  did.  [Ref.  35]  The 
inverse  ARMA  program  is  in  Appendix  H. 

2.  Performance 

The  effectiveness  of  the  inverse  ARMA  was  tested  against  the  synthetic 
scatterer  model  discussedsin  Chapter  II.  This  is  the  only  available  signal  for  which 
the  feed-forward  coefficients,  bf^  were  known.  For  actual  signals,  the  CadzOw- 
Solomon  algorithm  would  have  to  be  modified  to  extract.these  coefficients.  Presently 
the  algorithm,  particularly  the  bias  compensation  portion,  is  optimized  to  deterniine 
the  %  coefficients;  the  feedback  coefficients  which  describe  the  target  resonances. 
The  requirement  to  know:both  the  feed-forward  and  feedback  coefficients  forces  the 
inverse  ARMA  identification  process  to  be  aspect,  dependent. 

The  results  of  the  an  inverse  ARMA  filter  processing  correctly  matched, 
as  well  as  five  percent  above  and  below,  synthetic  signals  is  summarized  in  Table  10. 
In  order  to  allow  a  comparison  of  these  results  to  those  of  the  K-Pulse  filter  acting 
on  these  same  targets,  the  late-time  integration  periods  were  kept  the  same  (5.14  ns 
for  the  high  Q,  1.64  ns  for  medium  and  low  Q).  Importantly,  though,  this  period  was 
slid  forward  in  time  to  approximately  0.742  nanoseconds  into  the  record,  well  into 
early-time.  Comparison  of  Table  10  with  Table  6  indicates  the  advantage  the  inverse 
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ARMA  process  realizes  in  discriminatjgn  using^early-time  ener^.  Figure  42  shows 
the  inverse  ARMA  filter  output  for  the  20  dB  medium  Q  case  and  Figure  43  shows 
the  output  for  the  five  percent  above  filter  output.  These  results  can  be  compared 
to  Figures  27  and  28  for  the  K-Pulse  filter.  Notice  that  the  correct  filter  output  is 
simply  a  reproduction  of  the  incident  pulse;  all  the  target  return;  has  been 
deconvolved  out  of  the  data.  Both  the  high  and  medium  Q  cases  show  the  inverse 
ARMA  with  significantly  greater  ratios  of  incorrect-to-correct  filter  late-time  energies 
than  the  normal  K-Pulse  filter.  In  fact,  the  inverse  ARMA  performance  would 
probably  be  better  if  the  integration  time  were  not  limited  to  the  same  length  as  the 
K-Pulse  filter  for  direct  comparison. 


TABLE  10.  RATIOS  OF  OUTPUT  ENERGIES  FOR  lARMA  FILTER 


Test  Signal 

Signal  to 
Filter  Match 

SNR 

CO 

20  dB 

15  dB 

10  dB 

7dB 

High  o' 

5  %  above 

12880 

57.70 

18.90 

6.63 

3.80 

5  %  below 

14310 

66.12 

21.07 

7.92 

4.58 

Medium 

5  %  above 

4420 

58.12 

18.75 

6.42 

3.63 

5  %  below 

4630 

65.22 

22.10 

8.10 

4.76 

LowQ^ 

5  %  above 

417 

6.08 

Z43 

135 

1.13 

5  %  below 

422 

7.56 

3.25 

131 

1.45 

Noic  1:  Integration  time  =  5.14  ns 
Note  2:  Integration  lim&=  1.64  ns 
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Figure  42.  Correct  lARMA  filter  output  at  20  dB  SNR 
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Figure  43,  Five  percent  above  lARNIA  filter  output  at  20  dB  SNR 
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3.  Implementation  Considerations 

The  most  significant  change  infusing  an  inverse  ARMA  filter  process  for 
target  identification,  over  the  use  of  a  K-Pulse.filter  system,  is  theTact  that  it  is  /zor 
an  aspect  independent  process.  Thus  the  feed  forward  coefficients  of  the  target 
transfer  function,  the  bfjs,  would  have  to  be  determined  for  a  family  of  target  aspects 
oPmost  concern  for  each  target.  The  RAF  concept  described  earlier  would  be  much 
more  complex,  as  the  target  aspect  would  have  to  be  determined  by  another  radar, 
then  the  response  would  have  to  be  routed  to  a  filter  for  that  aspect  for  each  possible 
target.  Additionally,  the  number  of  filters  required  to  cover  the  target  aspects  of 
interest  (say,  a  60  degree  sector  about  the  nose)  would  have  to  be  determined.during 
the  extraction  process.  Perhaps  the  coefficients  for  certain  discrete  aspects  within  this 
region  could  be  determined  and  some  function  fit  to  them  to  determine  those  in. 
between.  In  any  case,  it  is  unlikely  that  the  feed-forward  portion  of  the  target 
transfer  function  would  be  exactly  matched  in  the  lARMA  filter.  This  would  then 
require  a  closed  loop,  adaptive  system  of  filters  that  could  start  with  the  coefficients 
as  extracted  in  the  filters  nearest  the  target  aspect  detected.  Then,  the  adaptive 
controller  would  dither  or  optimize  the  filter  coefficients  in  real-time  to  draw  out  the 
greatest  ratio  of  energies  output  from  the  filters.  This  would  require  much  increased 
integration  time  on  the  target  and  an  extremely  fast  and,  perhaps,  expert  adaptive 
controller.  [Ref.  35]; 

The  concerns  expressed  above  for  implementation  of  the  K-Pulse  filter, 
such  as  frequency,  bandwidth,  and  signal-to-npise  ratio,  remain  for  the  inverse  ARMA 
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concept.  Considering  the  increased  complexity  due  to  the  dependence  on  target 
aspect,  practical  implementation  of  the  inverse  ARMA  filter  process  for  target 
identification  would  significantly  lag  that  of  a  FIR-based,  aspe^ct  independent  method 
such  as  the  K-Pulse  filter. 
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V.  SUMMARY  AND  CONCLUSIONS 


A.  SUMMARY 

This  thesis  has  examined  a  resonance-based,  aspect  independent  radar  target 
identification  system  to  include  both  major  subsystems:  pole  extraction  and 
annihilation  filtering.  Chapter  I  presentedithe  theoretical  analysis  of  electromagnetic 
scattering  which  gives  rise  to  target  resonances  based  solely  on  structure  and 
composition,  independent  of  incident  field  aspect  or  polarization.  Chapter  11 
discussed  the  four  classes  of  signals  used  to  evaluate  the  system’s  algorithms,  namely: 
ARMA  generated  synthetic  data,  TDIE  thin-wire  data,  actual  measured  thin-wire 
data,  and  data  collected  from  scale  modehaircraft.  Chapter  III  presentedrthe  theory 
of  the  leading  candidate  pole- extraction  algorithm:'the  Cadzpw-Solomon  algorithm, 
results  of  its  use  on  the  target  signals,  and  some  initial  implications  for  real  world 
application.  Chapter  IV  applied  those  extraction  results  in  resonance  annihilation 
filters  to  evaluate  target  discrimination  and  discussed  relevant  areas  of  concern  for 
actual  implementation  of  this  portion  of  the  system.  Additionally,  the  inverse 
ARMA,  an  aspect  dependent  resonance-based  target  identification  scheme  was 
compared  to  the  aspect  independent  method. 


B.  CONCLUSIONS 


The  eventual  goal  of  application  of  a  resonance-based  target  identification 
system  is  to  gain  early  identification  of  non-cooperative  aircraft  targets.  The 
scattering  of  these  class  of  targets,  as  modeled  by  the  scale  aircraft  used  in  this  thesis, 
is  extremely  complex  and  spectrally  rich,  as  compared  with-that  of  the  more  simple 
thin-wire  targets  used  herein.  Intermediate  levels  of  target  complexity,,  such  as  few 
element  thin-wire  structures  should  be  investigated  via  these  algorithms  in  order  to 
refine  and  mature  their  application. 

The  Gadzow-Solomon  pole  extraction  algorithm  does  adequately  extract  poles 
of  simple,  thin  wire  targets,  both  simulated  and  actual.  It  requires  a  signal-to-hoise 
ratio  above  20  dB  and  its  result  is  sensitive  tb  its  various  input  parameters,  especially 
the  starting  point  of  its  input  data  relative  to  the  end  of  early-time.  The  over¬ 
estimation  of  the  system  order  and  resultant  ill-conditioning  of  the  data  matrix,  while 
necessary  to  allow  for  noise,;does  introduce  some  bias  in  the  extracted  poleiocations. 
The  complexity  of  the  frequency  spectrum  of  the  scattering  from  more  complex 
targets,  such  as  the  model  aircraft,  serves  as  a  limiting  factor  on  the  algorithm’s  ability 
to  extract  poles. 

Implementation  of  a  pole  extraction  scattering  facility  is  feasible  with  the  major 
caveat  being  accessibility  to  relevant  targets  for  scattering  analysis.  The  choice 
presented  is  between  actual  full-scale  targets,  which  will  be  difficult  to  obtain  and 
work  with,  considering  their  size;  and  scale  model  targets,  whose  resonance  frequency 
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range  is  more  amenable,  but  whose  fidelity  to  the  full-scale  targets  must  be 
thoroughly  evaluatedibefpre  implementation. 

The  K-Pulse  resonance  annihilation  filter  does  effectively  discriminate  simple 
targets  with  sets  of  poles  separated  by  five  percent  in  the  presence  of  noise  down  to 
10  dB  SNR.  However,  more  complicated  targets  present  an  increasingly  difficult  task 
because  their  pole  locations  are  many  and  are  closely  arid  irregularly  spaced.  The 
K-Pulse  filter,  in  its  simplest  foriri,  has  a  finite  limit  on  the  frequency  resolution  of 
zeros  of  different  targets  it  seeks  to  discriminate.  This  is  due  to  the  finite  width  of 
a  zero  null  in  the  frequency  response  of  the  filter. 

Implementation  of  the  resonarice  annihilation  filter  radar  for  fighter-sized 
targets  will  require  operation  in  the  medium  HF  to  low  VHF  bands:  approximately 
7  to  100  MHz.  The  design  challenges  here  are  the  same  as  for  ariy  radar  system  in 
this  regime  and  with  the  additional  requirement  to  obtain  as  precise  a  spatial  arid 
temporal  resolution  as  possible  in  order  to  isolate  ;the  target  return  from  noise, 
clutter,  and  other  targets.  Analysis  of  the  range  equation  for  a  system  operatingjat 
300  nautical  miles  indicates  the  requirement  for  250  kW  of  average  transmit  power. 
These  requirements  are  within  the  state-of-the-art  for  radar  systems.  The  greatest 
technological  emphasis  required  will  be  in  the  application  of  advanced  coniputingcor 
artificial  intelligence  to  control  the  radar  system  precisely  and  assist  the  operator  in 
analysis  of  the  returns. 

The  inverse  ARMA  algorithm  is  superior  to  the  K-Pulse  algorithm,  but  is  more 
difficult  to  implement  in  that  it  is  not  aspect  independent.  This  would  require  a 
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modification  of  the  pole  extraction  process  to  allpv.  extraction  of  poles  and  zeros  of 
the  target’s  transfer  function.  This  aspect  dependence  also^addsJthe  requirement  to 
create  many  filters  for  a  single  target  to  account  for  aspect-  This  additional 
requirement  would  further  emphasize  high  speed  processing  to  allow  tailoring  o&the 
filters  to  account  for  the  feed  forward,  aspect  and  excitation-depehdent^portion  of.the 
target  response. 

C.  AREAS  FOR  FURTHER  RESEARCH 

The  identification  of  non-cooperative  aircraft  targets  may  be  accomplished  by 
other  means,  such  as  high-range  resolution  or  synthetic  aperture  radars,  or  other 
target  characteristics.  As  all  these  technologies  mature,  the  relevant  benefits  and  costs 
of  each  must  be  compared. 

The  specific-algorithms  discussed  earlier  for  resonance^based  identification  must 
be  refined  by  application  to  more  intermediate  levels  of  target  complexity. 
Additionally,  study  of  the  frequency  response  of  these  less  complicated  targets,  as  well 
as  aircraft  targets,  may  discern  the  appropriate  bandwidth  requirements  to  allow  the 
most  effective  target  discrimination  possible. 

D.  UTILITY 

As  the  technology  of  the  attacker  becomes  more  advanced,  effective,  and 
prolific,  i.e.  stealth,  longer  range  precision  guided  munitions,  etc.,  the  technology  of 
NCTR  becomes  more  imperative.  Resonance-based  target  identification  is  such  a 
technology  that  holds  promise  for  development  into  a  viable  means  of  NCTR. 
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APPENDIX  A.  A^IA  COEFFICIENT  GENERATOR 


A.  PROGRAM  DESCRIPnON 

This  program  generates  the  and  b,.  coefficients  of  a  scatter  er  transfer  function 
given  by 


mz)  = 


Bq  +  b^z~^  +  ♦  •  •  +  bjZ' 
1  +  +  •  *  •  + 


after  reading  a  data  file  containing  the  ^-plane  poles  for  each  of  the  pole  :pairs. 
The  program  also  reads  in  the  sampling  inforrr  ation  to  allow  conversion  of  the  input 
j-plane  values  to  the  z-plane.  Additionally,  the  program  allows anput  of  a  complex 
multiple  of  each  z-plane  pole  pair  to  allow  simulation  of  the  relative  amplitude  and 
phase  of  the  transfer  function  poles. 

This  program  was  written  by  Capt.  Ti  J.  Murphy  and  Capt.  P.  C.  Reddy  arid 
uses  the  subroutine  POLY  written  by  Prof.  M.  A.  Morgan. 


PROGRAM  LISTING 


PROGRAM  HCOEF 

REAL*8  T,  A(0:100),  B(0:2),NUM1  (0:100).ATMP.'0:1 00),NUM2(0;1 00) 
REAL*8  SPOLR(30)i  SPOLI(30),Pl,MAG(30)iPHASE(30),C{30),D(30) 
INTEGER  NPTS,  NR,  N 
COMPLEX*=16  SPdL(30),  ZPCjL(30), 

CHARACTER*16  FNAME,  PFNAME 

Rl=3.141 59265 

WRlTE(*,*) ’Enter  fileriam.7  for  s-plane  poles  ’ 

READ{*,120)  PFNAME 


CPEN(1,FILE=PFNAME); 

READ(1,130)  NR 
WRrTE(*,*)  ’Enter  time  interval’ 

READ(*,*)  T 

WRITE(*,*)  ’Enter  number  of  points’ 

READ(^,*)  NPTS 
DO  10  1=  1,NR 
R^D(1  ,*)  SPOLR(I) 

READ(1,*)  SPOU(I) 

SPOL{l)=DCMPLX(SPOLR(l),SPOLI{l)) 

ZPOL(I) =CDEXP(DCMPLXCr/FLOAT(NPTl^  .))*3P01;(I)) 
wr|te(*,*)  ’Zrplane  pole  =’,2pol(0 
WRITE(*,'*)  ’Enter  magnitude  of  pole’,! 

READf*,*)  MAG(I) 

WHITE(*,*}  ’Enter  phase;(in  degrees)  of 
READ(^ ,*)  PHASE(l) 

PHASE(l)=0.0 
MAG(I)=1.0 

C  CONVERT  MAG  AND  PHASE  TO  RECTANGULAR  OOORD 
C(l)=MAG(i)*COS(PHASE(l)) 

D(I)=MAG(I)*SIN(PHASE(I)) 

10  CONTINUE 
A(0)=1.0 

A(1)=-2.0*REAL(ZPOL(1)) 

A(2)=(REAL(ZPOL(1)))**2+(AIMAG(ZPOL(1)))**2 

NUM1(0)=2.0*C(1) 

NUM1(1)=-2.0*(REAL(ZPOL(1))*C(1)+AIMAG(ZPOL(1))*DC?.)) 

NUM1(2)=0,0 

N=2 

DO  15  1=2, NR 
NUM2(0)=2.0*C(1) 

NUM2(1)=-2.0*(REAL(ZPOL(I))*C(I)+AIMAG(ZPOL(I))*D(I)) 

NUM2(2)=0.0 

B(0)=1.0 

B(1)=-2.0*REAL(ZPCL(I)) 

B(2) = (REAL(ZPOL(i)))**2+ (AIMAG{ZPOL(l)))**2 

DO  17K=0.N 

ATMP{K)=A(K) 

17  CONTINUE 

0  CALCULATE  NEW  DENOMINATOR 
CALL  POLY(A,B,N) 

C  CALCULATE  NEW  NUMERATOR 
N=N-2 

CALL  POLY(NUM1,B,N) 

N=N-2 

GALL  POLY(ATMP,NUM2,N) 

DO  18  L=0,N 

NUM1  {L)=NUM1  (L)+ATMP(L) 

18  CONTINUE 
15  CONTINUE 
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WRITE(*;^)  ’Enter  name  dWenominator  coefficient  outputsfile’ 

READ(*r120)^NAME 

OPEN(3;FlLE==t-NAME) 

WRrrE(3,130)  N 
DO  20  1=1, N 
WRrrE(3,*)  A{l) 

20  CONTINUE 

CL6SE(3) 

WRrTE(^,*)  ’Errter  name  of  numerator  coefficient  oulpui  file’ 

READ('‘,126)  FNAME 

OPEN(3iFlLE=FNAME) 

WRrrE{3,130)  N 
1^0  21  1=0, N 
VVRITE(3,*)  iMUMI(i) 

21  CONTINUE 
CLOSE(3) 

110  FORiVIAT(E16.10) 

120  FGRMAT(A) 

130  FORMAT(l5) 

END 

C 

SUBROUTINE  POLV(A,B,N) 

C 

C  Multiplying  {B(0)2**2  +,B(1)*z  +  B(2)}  x 
C  {A(0)z**N  +  A(1)*z**(N-1)  +  A(2)*2**{N-2)  +  ....  +  A(N)}  = 

C  C(0)z**(N+2)  +  C(1)*z**(N+1)  +  C(2)*2**N  -h  ....  +  C(N+2) 

C 

C  -Computing  C(n)  co#icients  and  storing  in  A(n)  while 
C  iricremehting  N  -->  N  +  2 
C 

REAL*8^A(0:1 00),B(G:2).C(0:1 00) 

C  A(0)=1.p 

C  Initialize  on  first  call  to  routine 
C  IF(N.GE;2)  GO  TO  11 

■C  N=2 

C  A(1)=E(1) 

C  A(2)=B(2) 

C  GO  TO  44 
11  G(0)=B(0)*A(0) 

C(1)=B(0)*A(1)+B(1)*A(0) 

DO  22  1=2, N 

22  C(I)=B(0)*A(I)+A(I-1)*B(1)+A(I-2)*B(2) 
C(N+1)=A(N)*B(1)+A(N-1)*B(2) 

C(N+2)=A(M)*B(2) 

N=N+2 
DO  33  l=0,N 
33  A(I)=C(I) 

44  CONTINUE 
RCTURN 
END 
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APPENDIX  B.  ARMA  SIGNAL  GENERATOR 


A.  PROGRAM  DESCRIPTION 

This  programs  generates  the  tu-ie  response  of  an  ARMA  system  via  the 
equation 


k=N  k=L 

3<”)  =  E  o^yin-k)  s. 

k=l  k=0 

where  x(«)  is  the  input  excitation  tecord,  a/,  and  b/.  are  the  coefficients  determined 
by  the  program  in  Appendix  A,  a>‘;d^y(/z)  is.the  output  data  record. 

This  program  was  written  by  Prof,  M.  A.  Morgan. 


B.  PROGRAM  LISTING 


Program  ARMA2 
C 

C  Computing  y[n]  response  for  N-th  order  ARMA  filter  due  to  x[n) 
C  input,  with  coefficients: 

C 

C  HR  a(k)fork=1,N 

C  FIR  b(k)  fork=0.L 

C 

C  by  M.A.  Morgan  5/23/90 

C 

REAL*8  x(0:2047),y(0:2047),a(1 :30),b(0:32) 

CHARACTER*64  TITLE, TITL 
characters  6  FNAME 
C  Entering  D.E.  Coefficients 

WRITE(*,*)  ’Enter  Filename  For  Recursive  a(k):’ 

WRITE(*,*)  ’Use  0  for  FIR  Filter' 

READ(*,1 00)  FNAME 
!F(FNAME.EQ.’0’)  GO  TO  15 
OPEN(2,FILE=FNAME) 

C  Assuming  that  a(iv)  are  reversed  polarity  in  file. 

READ(2,1i6)  N 
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write(*,^)  ’Write  k,  a(k):’ 

DO  i1ik=1,N 
READ(2j*)  a(k) 
a(k)=-a(k); 
write(*j*):k,a(k) 

11  CONTINUE 
CLOSE(2) 

15  CONTINUE 

WRITE(*;^)  ’Enter  Filename  for’Non-Recursive  b(k):’ 

READ(*,100)  FNAME 

OPEN(2,FrLE=FNAME) 

READ(2,T10)  L 
write(^,*)  ’Write  k,  b(k):’ 

DO22-k=0,L 
READ(2*);b(k) 
write(*,*)  k,b(k) 

22  CONTINUE 
CLOSE(2) 

WRITE(**)  ’Enter  Filename  for  x[n]  Plot  File:’ 
READ(*,1  00)  FNAME 
0PEN(2,FILE=FNAME) 

R^D(2,100)  TITL 
READ(2.110)  NX 
READ(2,120)'XQ 
READ(2,120)  XQ 
wfite(*,*)  ’Write  m,  x(m):’ 

D0  24m=0,NX-1 
READ(2,120)-x(m) 
write(^,*).  m,x(m) 

24  CONTINUE 
CL0SE(2) 

C  Setting  up  Plot  File  for  y(n] 

WRrTE(*,*)  ’Enter  Number  of  Points:  (.le.  2048)’ 
READ(**)  NY 

WRITE(*,*)  ’Enter  time  interval  (nanoseconds)’ 

REAp(**)Tmax 

Tmin=0.0 

WRITE(*,*)  ’Enter  filename  for  y[n]  response  output’ 
READ(MOO)  FNAME 
TITLE=  ’ARMA  Filter  y[n]  Response’ 
OPEN(1iFILE=FNAME) 

WRITE(1,100)  TITLE 
WRITE(1,110)  ,NY 
WRITE(1,i20)  Tmin 
WRITE(1,120)  Tmax 

C  Initializing  then  iterating  to  form  y{m] 
y(0)=b(0)*x(0) 
write(1,120)  y(0) 

D0  44m=1,NY-1 
y(m)=0.0 
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Lmax=MIN(m,L) 

DO  27  k=0,Lmax 
27  y(m)=y(m)+b(k)*x(m-k) 

Nmax=MIN(m,N) 

DO  33  k=1,Nmax 
33  y(m)=y(m)+a(k)*y{m-k) 

wiite(*,*)  y(m) 

WRITER,  120)  y(m) 

44  CONTINUE 

WRITE(*,*)  Tbats  All  Folks  !‘ 
CUOSEO) 

100  FORMAT(A) 

110  F0RMAT(I5) 

120  FORMAt(E12.6) 


APPENDIX  C.  NOISE  POLLUTION  GENERATOR 

A.  PROGRAM  DESCRIPTION 

This  program  reads  in  a  data  file  and  computes  the  signal  power.  After  inputing 
a  seed  between  zero  and  one,  it  uses  a  routine  from  the  FORTRAN  Scientific 
Librafay  called  STRNUM  to  generate  numbers  under  a  zero  mean  uniform 
distribution.  This  routine  is  called  1000  times  for  each  additive  noise  point 
corresponding  to  each  point  of  the  input  record,  in  order  to  drive  the  random 
number  to  a  zerormean  Gaussian  distribution.  The  variance  of  the  fandom  noise 
generated  is  scaled  to  accomodate  a  user  input  signal-to-noise  ratio. 

This  program  is  based  oh  a  sipal  generator  originally  written  13y  Prof..  M.  A. 
Morgan,  modified  succe.ssively  by  Lt.  L.  A.  GheekSj  Lt.  M.  S.  Simon,  and  Capt.  T.  J. 
Murphy. 


B.  PROGRAM  LISTING 


PROGRAM  NPLUS.FOR 
C 

C  PROGRAM  NOTES: 

C  1)  GAUSSIAN  NOISE  ADDITION  TO  SYN  SIGNAL  OPTION 
C  2)  PROGRAM  REQUIRES  8087  MATH  CO-PROCESSOR  TO  RUN 
C  3)  ORIGINAL  SYNGEN  PROGRAM  WRiTTEN  IN  BASIC 
C  BY  DR.  M.  A.  MORGAN 

C  4)  STW087P0R  PROGRAM  WRITTEN  IN  MSFORTRAN  77 
C  BY  L  A.  CHEEKS 

C  5)  UPDATE  AUG89.BY  LT  M.S.  SIMON 
C  6)  UPDATE  MAR90  BY  CAPT  T.  J;  MURPHY 
C 

C  ‘‘‘DECLARATIONS*** 
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c 

REAL  T9.T0,A;C9,PI.X0SUM,GSUM,GSUMPRM,R,CDB,X0(1 024), R0(100)jl0(10 
+0).A0(1 00)iPd(1 00),G(1 024), X1  (1 024),R,TZ/0.0/,DT1  ,DT,L2,X02SyM,X02 
+(1024),GADJ,GMEAN 
INTEGER  N,Z,I9,I1,NC0UNT 
CHARACTER;FN*1 6.TITLE*1 6,  ANS*1  ,HEADER*64 
DATA  PI/3.1415927/ 

C 

C  ***MAIN  BODY  OF  PROGRAM*** 

C 

C  ***DATA  INPUT  ROUTINE*** 

0 

170  WRITE(*,*)  'ENTER  NAME  OF  DATA  FILE  TO  BE  CORRUPTED' 

READ  (*,280)  TITLE 
OPEN(1,FILE=TITLE) 

READ  (1,260)  HEADER 
READ  (1,240)  NPTS 
READ=(1,250)  TBEGIN 
READ  (1,250)  TEND 
C 

DO  190  1=1, NPTS 
READ  (1,250)^X0(I) 

190  CONTINUE 
CLOSE(I) 

C 

WRITE(*,*)'DONE  READING  FILE' 

C 

X0SUM=0.0 

DO  80  1=1, NPTS 

X0SUM=X0SUM+(X0(I)**2) 

80  CONTINUE 
C 

C  XOSUM/NPTS  IS  AVERAGE  POWER 
X0SUM=X0SUM/NPTS 
C 

C  ***ROUTINE  TO  ADD  GAUSSIAN  NOISE  TO  SIGNAL*** 

C 

GSUM=0.0 

GADJ=0.0 

GMEAN=0.0 

120  WRITE(*,*)'ENTER  SEED  FOR  RANDOM  NUMBER  GENERATION  FOR  NOISE' 
WRITE(*,*)'(0  <  R  <  1.0)  0.21289  GENERATES  A  GOOD  DISTRIBUTION' 
READ(*,250)  R 
DO  130  1=1, NPTS 
CALL  CLMT(R,P) 

G(I)=P 

GADJ=GADJ+G(I) 

130  CONTINUE 

C  ADJUST  DISTRIBUTION  SO  MEAN  IS  0 
GMEAN=GADJ/FLOAT(NPTS) 
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DO  135  l=1,NPTS 
G(I)=G(I)-GMEAN 
GSUM=GSUM+(G(I)**2) 

135  CONTINUE 

C  GSUM/NPTS  IS  WERAGE  POWER  OF  GAUSSIAN  SIGNAL 

GSUM=GSUM/NPTS 

WRITE{**):?ENTER  DESIRED  SIGNAL  TO^NOISE  RATIO  IN  DB’ 
READ(**)  DB 

GSUMRRM=XOSUM/(10.0**(DB/10.0)) 

CDB=10.0*ALdG10(X0SUM/GSUMPRM) 

WRITE(* ,*)  'CHECK  DB  =  '.CDB 
RATIO=SQRT(GSUMPRM/GSUM) 

DO  140I=1,NPTS 
G(l)=RATIO*G(l) 

140  CONTINUE 

DO  150I=1,NPTS 
X0(l)=X0(l)+G(r) 

150  CONTINUE 
C 

C  ***ROUTINE  TO  NORMALIZE  SYN  GEN  SIGNAL*** 

C 

155  X02SUM=0.0 

DO  156I=1,NPTS 

X02(I)=X0(I)**2 

X02SUM=X02SUM+X02(I) 

156  CONTINUE 
L2=SQRT(X02SUM) 

DO  16bl=1,NPTS 
X0(I)=X0(I)/L2 

160  CONTINUE 

C 

C 

C  ***ROUTINE  TO  STORE  SYN  GEN  SIGNAL  DATA  IN  DATA  FILE*** 


WRITE(*,*). 'ENTER  NAME  OF  DATA  FILE  TO  CONTAIN  CORRUPTED  SIGNAL' 
READ  (*,280)  FN 

WRITE(*,*)>'ST6RING  corrupted  signal  data  in  data  FILE' 
OPEN(1,FILE=FN) 

WRITE(1,260)  HEADER 
WRITE(1,240)  NPTS 
WRIT^(1,250)  TBEGIN 
WRITE(1 ,250)  TEND 
C 

DO  17  1=1, NPTS 
WRITE(1,250)  X0(l) 

17  CONTINUE 
CLOSE(I) 

C 

240  FORMAT(l5) 
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250  FORMAT(E12.6) 

260  FORMAT(A64) 

270  FORMAT(A8) 

280  FORMAT(A16) 

WRITE(**)  ’CHECK  NEW  FILE(S)’ 

GO  TO  170 
STOP 
END 
C 

C  ***SUBROUTINES*** 

SUBROUTINE  CLMT(R,P) 

C 

C  THIS  SUBROUTINE  IS  BASED  ON  THE  CENTRAL  LIMIT  THEOREM  AND 
C  GENERATES  A  SEQUENCE  OF  INDEPENDANT  NUMBERS  WHICH  ARE 
C  NORMALLY  DISTRIBUTED  ABOUT  THE  INTERVAL  0  TO  1 
C 

REAL  P.R 
C 

P=0.0 

DO  170  1=1,1000 
CALL  STRNUM(R) 

P=P+R 

170  CONTINUE 

p=p/iooo;o 

RETURN 

END 

C 

SUBROUTINE  STRNUM(R) 

C 

C  THIS  SUBROUTINE  GENERATES  A  SEQUENCE  OF  NUMBERS  WHICH^ARE 
C  RANDOMLY  AND  UNIFORMLY  DISTRIBUTED  OVER  THE  UNIT  INTERVAL 
C  REF;  P.  280  OF  FORTRAN:  Scientific  Subroutine  Library 
C 

C  NOTES: 

C  1)  R  =  0.21289  GENERATES  A  GOOD  UNIFORM  DISTRIBUTION 
C 

BB=1.0 

P1=R*317.0 

R=AMOD(P1,BB) 

■RETURN 

END 


APPENDIX  D.  K-PULSE  GENERATOR 


A.  PROGRAM  DESCRIPTION 

This  program  creates  an  K-Pulse  filter  impulse  response  after  reading  a  file 
containing  the  j-plane  poles  to  be  annihilated.  This  program  allows  scaling  of  the 
filter  impulse  response  as  welhas  integer  delays  bewteen  filter  weights,  each  to  allow 
additionahdegrees  of  freedom  in  the  filter  design.  In  this  thesis,  these  parameters 
were  set  to  one. 

This  program  was  written  by  Prof.  M.  A.  Morgan,  with  susequent  modifications 
by  Lt.  M.  S.  Simon  and  Capt.  P.  C.  Reddy. 


B.  PROGRAM  LISTING 


PROGRAM  :KP 
C 

C  PROGRAM  NOTES 
C  1)  program' GENERATES  KILL  PULSES 

C  2)  KEYBOARD  INPUTS:  POLES/REAL, IMAG 

0  3)  OUTPUT  TO  DATA  FILE:  K-PULSE  TIME  SERIES 

C  4)  PROGRAM  REQUIRES  8087  MATH  CO-PROCESSOR 
C  5)  PROGRAM  WRITTEN  IN  MSFORfRANBY  PROF.  M.A.  MORGAN 
C  6); INTERACTIVE  MODS  BY  LT  M.S.  SIMON 

C  7)  ADDITONAL  MOS  BY  GAPT  P.C.  REDDY 
C 

C  K-Pulse  Generation  Using  Z-transform  product  Non-Recursive 
C  Filter  Design  Approach  by  M.A.  Morgan  8/24/89. 

C 

REALSR(20),SI(20),A(0:100),B(2),H(401),SF 
CHARACTER  FNAME*1 6,TITLE*64,ANS*1  ,PFNAME*1 6 
WRlTE(*, *)  ’***  K-PULSE  GENERATOR  ***’ 

WRITE(*,*) 

1  WRITE(*,*): ’ENTER  NUMBER  OF  POLE  PAIRS  (<=50):’ 

WRITE(*,*)=!ENTER  SCALE  FACTOR' 

READ{**)SF 


WRITE(* ,*)  ’ENTER  TIME  WINDOW  (NANOSEC);’ 

READ(* ,*)  T0= 

WRrrE(*  *)  ’ENTER  NUMBER  OF  TIME  POINTS:’ 

READ(*,*)  NT 
DT=T0/(NT-1;0) 

WRrrE(V)  ’ENTER  TIME-STEPS  BETWEEN  FILTER  WEIGHTS  (>=1):’ 
READ(* ,*)  KS 

WRITE(*i*)  ’ENTER  FILENAME  OF  SPLANE  POLES’ 

READ(*,100)  PFNAME 
OPEN(2,FILE=PFNAME) 

READ(2 ,*)  NP 
DO  11  N=1,NP 
READ(2,*)  SR(N) 

READ(2.*)  SI(N) 

11  CONTINUE 
CLOSE(2) 

0 

C  INPUT  DATA  CHECK 
C 

WRITE(*,*)  ’  ’ 

WRITE(*,?)’INPUT  DATA  CHECK’ 

WRITE(*  *)’NUMBER  OF  POLE  PAIRS  IS  ’,NP 
WRITE(*,*)’SGALE  FACTOR  IS  ’.SF 
WRITEh*)’TIME  WINDOW  IS  ’,T0,’  NANOSEC’ 

WRITE(*i*)’NUMBER  OF  TIME  POINTS  IS", NT 

WRITE(*,*)’NUMBER  OF  TIME  STEPS  BETWEEN  FILTER  WEIGHTS  IS  ’,KS 
WRITE(*,*)’  ’ 

DO  12N=1,NP 

WRITE(*,*)’SIGMA’,N,’  =  ’,SR(N) 

WRITE(*.*)’OMEGA  ’,N,’  =  ’,SI(N) 

PAUSE 

12  CONTINUE 

13  WRITE(*,*)’IS  DATA  CORRECT?’ 

READ(*,200)ANS 

IF((ANS.EQ.’Y’).OR.(ANS.EQ.’y’))THEN 
WRITE(*,*)’DATA  CORRECT 
GOTO  14 

ELSEIF((ANS.EQ.’N’).OR.(ANS.EQ.’n’))THEN 
WRITE(*,*)’DATA  ERRORS.  START  AGAIN’ 

GOTO  1 
ELSE 

WRITE(*,*)’TRY  AGAIN!  Y  OR  N’ 

GO  TO  13 
ENDIF 
C 

C  SCALING  POLES 
C 

14  D0  15N=1,NP 
SR(N)=SR(N)/SF 
SI(N)=SI(N)/SF 
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15  CONTINUE 
C: 

C  CALCULATIONS 
C 

MP=0 

DO  22  N=1,NP 
E=EXP(KS*SR(N)*DT) 

C=COS(kS*SI(N)*DT) 

B(1)=-2.0*E*C 

B(2)=E*E 

CALL  P0LY(A,B.MP) 

22  CONTINUE 
WRITE(**) 

WRITE(*;^)  ’UN-NORMALIZED  K-PULSE  FILTER  WEIGHTS:’ 
WRITE(**) 

D0.33  M=0,MP 
WRITE(* ,*)  M,A(M) 

33  CONTINUE 
PAUSE 
C 

C  NORMALIZE  K-PULSE  TO  UNIT  ENERGY 
C 

ENERGY=A(0)**2 
M1=KS*MP-H 
D0  44M=1,MP 

44  ENERGY=ENERGY-FA(M)**2 
RMS=SQRT(ENERGY) 

C 

C  DEFINING  UNIT  SAMPLE  RESPONSE 
C 

H(1)=A(0)/RMS 

1=1 

DO  55N=1,MP 
M=0 

50  M=M-t-1 

l=l-H 

IF(M.GE.KS)  GO  TO  55 

H(l)=0.0 

GO  TO  50 

55  H(I)=A(N)/RMS 
C 

C  STORING  K-PULSE  IN  DATA  FILE 
C 

TS=0.0 

TF=KS*MP*DT 

WRITE(*,*)  ’STORING  K-PULSE  IN  DATA  FILE’ 

WRITE(*,*) 

WRITER, *)  ’ENTER  OUTPUT  FILE  (<=  16  CHAR):’ 

READ(MOO)  FNAME 

WRITE(*,*)  ’ENTER  HEADER  FOR  OUTPUT  FILE  (<=  64  Char):’ 
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READ(M  00)  TITLE 
OPEN(1,FILE=FNAME) 

.WRITE(T,100).TITLE 
WRjtE(1,;i10)  M1 
WRltE(ni2(^;TS 
WRrrE(T,120)  TF 
D0  66N=1,M1 
WRITE(1,120):H(N) 

66  CONTINUE 

88  WRITE(*,*)’K-PULSE  GENERATOR  PROGRAM  COMPLETE’ 

WRITE(*  *)’CHECK  NEW  FILE’ 

100  FORMAT(A) 

110  F0RMAT(I5) 

120  FORMAT(E12.6) 

200  FORMAT(A1) 

STOP 

END 

C 

SUBROUTINE  POLY(A,B,N) 

C 

C  MULTIPLYING  {2**2  +  B(1)*z  +  B(2)}  x 
C  {2**N  +  A(1)*z**(N-1)  +  A(2)*2**(N-2)  +  ....  +  A(N)}  = 

C  z**(N+2)  +  C(1)*2**(N+1)  +  C(2)*2**N  +  ....  +  C(N+2) 

C 

C  COMPUTING  C(N)  COEFFICIENTS  AND  STORING  INTO  A(N)  WHILE 
C  INCREMENTING  N  ->  N  +  2 
C 

REAL  A{0:100),B(2),C(0:100) 

A(0)=1.0 

C  INITIALIZE  ON  FIRST  CALL  TO  ROUTINE 
IF(N.GE.2)  GOTO  11 
N=2 

A(1)=B(1) 

A(2)=B(2) 

GO  TO  44 

11  C(1)=A(1)+B(1) 

DO  22  1=2, N 

22  C(I)=A(I)+A(I-1)*B(1)+A(I-2)*B(2) 

C(N+1)=A(N)*B(1)+A(N-1)*B(2) 

G(N+2)=A(N)*B(2) 

N=N+2 
DO  33  1=1, N 
33  A(I)=C(I) 

44  CONTINUE 
RETURN 
END 
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APPENDIX  E.  DOUBLE  GAUSSIAN  SMOpTKING  FUNCnON  GENERATOR 

A.  PROGRAM  DESCRIPTION 

This  program  creates  a  double  Gaussianismoothihg  function  used  to  allow  the 
use  of  the  difference  equation  based:  K-Pulse  filter  in  the  presence  the  noise.  The 
program  uses  input  values  of  the  ten  percent  height  of  the  low  and  high  frequency 
ends  of  the  double  Gaussian  frequency  response  to  determine  the  pulse  widths  in  the 
time  domain.  The  user  defines  a  threshold  value  to  truncate  the  tails  of  the  double 
Gaussian  function  to  limit  its  temporal  extent.  The  .program  uses  the  PLTSUB 
subroutine  written  by  Prof.  M.  A.  Morgan,  which  uses  the  Grafmatics  FORTRAN 
graphics -package,  to  create^a  plot  of  the  result. 

This  program  was  originally  written  by  Prof.  M.  A.  Morgan  and  was  modified- 
by  Capt.  P.  C.  Reddy. 


B.  PROGRAM  LISTING 

PROGRAM  DGSMF 
C 

C  Adaptive  Double-Gaussian  Smoothing  Function  Generator 
C  By  M.A.  Morgan  2/22/90/=Modified  by  P.C.  Reddy  8/1/90 
C 

REAL  Y1  (1 024),DT,T1  ,T2,F1  ,F2,TMIN.TMAX,TMAX1  ,THR 

REAL*8  Y(1 024), El  ,E2,F(1 024),  A1  ,A2,ALFA1  ,ALFA2,C1  ,C2 

INTEGER  Npfs,  N 

CHARACTER*1  YN 

CHARACTER*64  FNAME,TITLE 

WRitE(*,*)  ’Enter  No.  Time  Points  ’ 

READ(*,*)  NPTS 

TMIN=0.0 

WRITE(*,*) 
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WRrTE(*,*)  ’Enter  Time  Window  in  NSgc:  ’ 

READ(*,*)  TMAX 
DT=TMAX/FLOAT(NPTS-i) 

44  CONTINUE 

WRiTE(*,*)  ’Enter  Upper  40%  Freq  Rolloff  in  GHz  ’ 

READ(*,«)  F1 

T1=1.47/F1 

WRITE(*.*)  ’Narrow-10%:Pulse  Width  in  NSec:=  ’,T1 
WRITE(*,*)  ’Enter  Lower  10%  Freq  RoHofFin  GHz  ’ 
READ(* ,*)  F2 
WRrrE(*,*) 

T2=.31/F2 

WRITE(*,*)  ’Wide  10%  Pulse  Width  in  NSec:=  ’,T2 
WRiTE(*,*) 

WRITE(*,*)  ’Are  These  Pulse  Widths  OK  ?  (Y/N)  ’ 
READ(M00)YN 

IF^N.EQ.’N’.dR.YN.EQ.’n’)  GO  TO  44 
WRITE(*,*)  ’Enter  THRESHOLD  ’ 

READ(*,*)  THR 
ALFA1  =4*DLOG(1 0)/T1  ’^*2 
ALFA2=4*DLOG(1  Q)tT2*^2 
A1  -  DSQRT(ALFA1  )/(DSQRT(ALFA1  )-DSQRT(ALFA2)) 
A2=A1-1.0 
1=0 

DO  10  K=1;NPTS 

T=(K-1.0-.(NPTS/2))*DT 

C1=ALFA1*T**2 

C2=ALFA2*t**2 

IF(C1.GT.THR)THEN 

El  =0.0 

ELSE 

E1=A1*EXP(-C1) 

ENDIF 

IF(C2.GT.THR)THEN 
F(K)=0.0 
GOTO  10 
ELSE 

E2=A2*EXP(-C2) 

ENDIF 

F(K)=E1-E2 

9  IF(F(K).NE.0.0)THEN 
1=1+1 

Y(I)=F(K) 

Y1(I)=SNGL(Y(I)) 

N=l 

ENDIF 

10  CONTINUE 
TMAX1=DPN 

C  EGA  Assumed: 

NS=1 
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C  Impact  Printer  Assumed: 

NFD=1 

C 

TITLE=’DOUBLE-GAUSSIAN  SMOOTHING  FUNCTION’ 
CALL  PLtSUB(TlTLE,N,TMIN,TMAX1  ,Y1  .NB.NFD) 
WRITE(*it)  ’Do  You  Want  to  Store  These  Results  ?'(Y/N) 
READ(MOb)YN 

IF(YN.EQ:’N’.OR.YN.EQ.’n’|GO  TO  99 
WRITE(*,*)  ’Enter  File  Narne  for  Plot  Data  Storage:  ’ 
READ(''.100)  FNAME 
OPEN(3,RLE=FNAME) 

WRITE(3.100)  TITLE 
WRITE(3,i;iO)  N 
WRITE(3,T20)  TMIN 
WRITE(3,120)  TMAX1 
DO  88  l=1iN 
WRITE(3,120)  Y(l) 

83  CONTINUE 
CL03E(3): 

99  CONTINUE 

WRITE(*,f')  ’Program  Completed’ 

100  FORMAT(A) 

110  F0RMAT(I5) 

120  FORMAT(E12.6) 

STOP 

END 


118 


ooo  o  ooooo 


APPENDIX  F.  SMOOTHED  K-PULSE  GENERATOR 


A.  PROGRAM  DESCRIPTION 

This  program  is  an  adapataion  of  a  linear  convolution  program  written  by  Prof. 
M.  A.  Morgan.  This  program  reads  two  data  files  and  convolves  them  and 
normalizes  the  output  waveforrn  energy  before  depositing  the;  result  in  user  specified 


data  file. 


B.  PROGRAM  LISTING 
PROGRAM  CVN 

Convolution  of  PLOT  Format  Data:Files  X(n).*Y(n)  to  Produce 
Z(n).  Provides  Local  Plotting  and  PLOT  Output  File, 

By  M.  Morgan,  8/21/89.  Modified  by  P.. Reddy,  8/25/90. 

REAL  X(1 024),  Y(1 024),Z(2048).7.1  f2048),TLATE 
CHARACTER*1  ANS,  BELL 
CHARACTER*8  FNX,FNY,FN2 
CHARACTER*64  TnX,TIPi',7lTZ 

BELL=CHAR(7) 

DATA  FILE  INPUT 

10  CONTINUE 

WRiTE(*,*)  THIS  PROGRAM  CONVOLVES  X{n)  *  Y(n)  =  Z(n)’ 
WRITE(*,*)  ’AND  PLACES  Z(n)  IN  USER  SPECIFIED  OUTPUT  FILE’ 
WRITE(*,*) 

WRITE(*,100)  BELL 

WRITE(*,*)  ’ENTER  NAME  OF  X(n)  INPUT  FILE’ 

V/RITE(*,*) 

READ(*,100)  FNX 
WRITE{*,*) 

WRITE(*,*)  ’X(n)  INPUT  FILE  IS  ’.FNX 
WRITE(*,*) 

WRITE(*,100)  BELL 

WRITE(*,*)  ’ENTER  NAME  OF  Y(n)  INPUT  FILE’ 


WRITE(*.*) 

READp,100)  FNY 
WRITE(**) 

WRlTE(*,?)  ■Y(h)  INPUT-FILE  IS  ’.FNY 
WRITE(*.*) 

C 

G  READING  INPUT  FILES 
C 

0PEN(1,FILE=FNX) 

0PEN(2,FILE=FNY) 

READ(1.100)  TITX 
READ(1,120)  NX 
READ(1,130)  rk\ 

READ(1,130)  TX2 

DO20;l=1.NX 

READ(1,130)X(I) 

20  CONTINUE 

WRITER,*)  ’X(t)  FILE  READ’ 

READ(2,100)  TITY 
READ(2,120)  NY 
READ(2.130)  TY1 
READ(2,130)  TY2 
DO  30  1=1, NY 
READ(2,130)  Y(l) 

30  CONTINUE 

WRITE(*,*)  ’YW  file  READ’ 

WRITE(*,*) 

CLOSE(I) 

CLOSE(2) 

C 

C  CONVOLUTION  ALGORITHM  V\/ITHOUT  ZERO  PADDING 
C  Written  by  M. A.  Morgan  5/19/89 
C 

WRITE(*,*)  ’CONVOLUTION  IN  PROGRESS’ 

WRITE(*,*) 

DO  50  N=1,NY 

Z(N)=0.0 

M2=MIN0(N,NX) 

DO  40M=1,M2 

Z(N)  =Z(N)  +X(M)*Y(N-M + 1 ) 

40  Z1(N)=Z(N) 

50  CONTINUE 
N1=NY+1 
N2=NX+NY-1 
DO  70  N=N1,N2 
Z(N)=0.0 
M1=N+1-NY 
M2=MIN0(N,NX) 

DO  60  M=M1,M2 
60  Z(N)=Z(N)+X(M)*Y(N-M+1) 
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70  GONTINUE 
C 

TZ1=TX1+TY1 
tZ2=TX2+TY2 
ENERGY=0.0 
DO  90:N=1,N2 

90  ENERGY=ENERGY+Z(N)*Z(N) 

DO  92  N=1,N2 

Z(N)  =Z(N)/SQRT(ENERGY) 

92  Z1(N)=Z(N) 

WRITE(*,*) 

PAUSE 

C 

C  ***  Plotting  Z(n)  *** 

•nTZ=’ConvoIution  Z(n)=X(n)*Y(n)’ 

CALL  PLTSUBCriTZ,N2,TZ1  ,TZ2,T71  ,TZ2,Z1 ,1 ,1 ,0) 

C 

WRiTE(MOO)  BELL 

WRITER,*)  ’ENTER  NAME  OF  Z(n)  FILE’ 

WRITE{**) 

READ(MOO)  FNZ 
WRITE(*,*) 

WRITE(**)  ’OUTPUT  FILE  IS  ’.FNZ 
OPEN(3.FILE=FNg 
WRITE(**)- 
WRITE(MOO)  BELL 

WRITE(*;^)  ’ENTER  HEADER  TITLE  FOR  Z(n)  FILE’ 

WRITE(**) 

READ(MOO)  TITZ 
WRITE(**) 

VVRITE(*  *)  ’SAVING  Z(t)  TO  OUTPUT  DATA  FILE’ 

WRITE(3.100)  TITZ 
WRITE(3.120)  N2 
WRITE(3,130)  TZ1 
WRITE(3,130)  TZ2 
DO  80i=1,N2 
WRITE(3,130)  Z(l) 

80  CONTINUE 
CLOSE(3) 

C  QUERY  TO  DO  ANOTHER  RUN 
C 

95  CONTINUE 
WRITE(**) 

WRITE(*,100)  BELL 

WRITE(*,*)  ’DO  YOU  WISH  TO  CONVOLVE  ADDITIONAL  WAVEFORMS?’ 
WRITE(**) 

READ(*,100)  ANS 

IF((ANS.EQ.’Y’).OR.(ANS.EQ.’y’))THEN 
GOTO  10 

ELSEIF((ANS.EQ.’N’).OR.(ANS.EQ.’n’))THEN 


121 


GO  TO  99 
ELSE 

WRITE(^;*)  TRY  AGAIN!  ■¥■  OR  'N'.’ 

GO  TO  95 
ENDIF 

99  CONTINUE 
WRiTE(*,100)iBELL 

WRiTE(* ,*)  ’CONVOLUTION  PROGRAM  IS  OVER.  CHECK  NEW  F!LE(S)’ 

100  FORMAT(A) 

120  FORMAt(l5) 

130  FORMAT(E12.3) 

STOP 

END 
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APPENDIX  G.  K-PULSE  CONVOLUTION  PROGRAM 


A.  PROGRAM  DESCRIPTION 

This  program  is  an  adapataion  of  a  linear  convolution  program  wriiten  by  Prof. 
M.  A.  Morgan.  This  program  reads  two  data  files  and  convolves  them  and  calculates 
the  average  energy  over  a  specified  time  of  the  output  record.  The  program  plots 
the  result  and  deposits  it  in  a  user  specified  data  file.  This  program  was  modified  b> 
Capt.  P.  C.  Reddy. 


B.  PROGRAM  LISTING 

PROGRAM  CV2 
C 

C  Convolution  of  PLOT  Format  Data  Files  X(n)*Y(n)  to  Produce 
C  Z(n).  Provides  Local  Plotting  and  PLOT  Output  File  Along 
C  With  Late-Time  Energy  Calculation. 

C  By  M.  Morgan,  8/21/89.  Modified  by  P.  Reddy,  8/25/90. 

C 

REAL  X(1024),Y(1024),Z(2048),Z1  (2048),TLATE 
CHARACTER*!  ANS,  BELL 
CHARACTER*8  FNX,FNY,FNZ 
CHARACTER*64  T1TX,TITY,TITZ 
C 

BELL=CHAR(7) 

C 

C  DATA  FILE  INPUT 
C 

10  CONTINUE 

WRITE(*,*)  THIS  PROGRAM  CONVOLVES  X{n)  *  Y(n)  =  Z(n)’ 
WRITE(*,*)  'AND  PLACES  Z(n)  IN  USER  SPECIFIED  OUTPUT  FILE’ 
WRitE(*,*) 

WRrrE(*,100)  BELL 

WRITE(*,*)  'ENTER  NAME  OF  X(n)  INPUT  FILE’ 

WRITEC,*) 

READ(*,100)  FNX 
WRITE(*.*) 

WRITE(*,*)  ’X(n)  INPUT  FILE  IS  ’.FNX 
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WRITE(**) 

WRITE(*,100)  BELL 

WRITE(*,*)  'ENTER  NAME  OF  Y(n)  INPUT  FILE’ 
WRITE(*.*) 

READ(MOO)  FNY 
WRITE(*,*) 

WRITE(* ,*)  ■Y(n)  INPUT  FILE  IS  '.FNY 
WRITE(*,*) 

C 

C  READING  INPUT  FILES 
C 

OPEN(1,FILE=FNX) 

OPEN(2,FrLE=FNY) 

READ{1,100)  TITX 
READ(1,120)  NX 
READ(1,130)  TX1 
READ(1,130)  TX2 
DO  20  1=1, NX 
READ(1,130)  X(l) 

20  CONTINUE 

WRITE(*,*)  ’XO)  FILE  READ' 

READ(2,100)  TITY 
READ(2,120)  NY 
READ(2,130)  TY1 
READ(2,130)  TY2 
DO  30  1=1, NY 
READ(2,130)  Y(l) 

30  CONTINUE 

WRITE(*,*)  ’Y(t)  FILE  READ’ 

WRITE(*,*) 

CLOSE(I) 

CLOSE(2) 

C 

C  CONVOLUTION  ALGORITHM  WITHOUT  ZERO  PADDING 
C  Written  by  M.A.  Morgan  5/19/89 
C 

WRITE(*,*)  'CONVOLUTION  IN  PROGRESS’ 

WRITE(*,*) 

DO  50N=1,NY 

Z(N)=0.0 

M2=MIN0(N,NX) 

DO  40  M=1,M2 
Z(N)=Z{N)+X(M)*Y(N-M+1) 

40  Z1(N)=Z(N) 

50  CONTINUE 
N1=NY+1 
N2=NX+NY-1 
D0  70N=N1,N2 
Z(N)=0.0 
M1=N+1-NY 
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M2=MIN0(N,NX) 

DO:60M=M1,M2 

Z(N) =2(N) +X(M)*Y(N-M + 1 ) 

60  Z1(N)=Z(N) 

70  CONTINUE 
C 

C  CALCULATING  LATE-TIME  ENERGY 
C 

83  TZ1=TX1+TY1 
TZ2=TX2+TY2 
WRITE(*.*) 

WRITE(**)  'ENTER  TL  (Nsec)  FOR  LATE-TIME  ENERGY  CALCULATION’ 
WRITE(**) 

READ(*  *)  TLATE 

WRITE(*,*): 'ENTER  TZ3  (Nsec)  FOR  END  OF’ 

WRltE(*.*)  'LATE-TIME  ENERGY  CALCULATION' 

WRITE(*,*) 

READ(*  *)  TZ3 

NLATE=NINT(TLATE*N2/CTZ2-TZ1))+1 
N3= NINT(TZ3*N2/(TZ2-TZ1 )) + 1 
IF(NLATE.LT.N2)  GO  TO  85 
WRITE(*,100)  BELL 

WRITE(*,*) 'TLATE  TOO  LARGE  -  TRY  AGAIN' 

GO  TO  83 
85  CONTINUE 

IF(NLA.TE.GE.1)  GO  TO  87 
WRITE(MOO)  BELL 

WRITE(*,*)  TLATE  TOO  SMALL  -  TRY  AGAIN’ 

GO  TO  83 
87  CONTINUE 
ENERGY=0.0 
DO  90  N=NLATE.N3 
90  ENERGY=ENERGY+Z(N)*Z(N) 

ENERGY=ENERGY/(N3-NLATE+1) 

WRITE(*,*) 

WRITE(*,*)  'LATE  TIME  ENERGY  =’,  ENERGY 
PAUSE 
C 

C  ***  Plotting  Z(n)  *** 

TITZ=’Convolution  Z(n)=X(n)*Y(n)’ 

CALL  PLTSUB(TITZ,N2,TZ1  ,TZ2,TZ1  ,TZ2,Z1 ,1 ,1 ,0) 

C 

WRITE(*,100)  BELL 

WRITE(*,*)  'ENTER  NAME  OF  Z(n)  FILE’ 

WRITE(**) 

READ(MOO)  FNZ 
WRrrE(*,*) 

WRITE(*,*)  'OUTPUT  FILE  IS  ’,FNZ 
OPEN(3,FILE=FNZ) 

WRITE(*,*) 
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WRITE(*,100)  BELL 

WRlTE(* ,*)  ’ENTER  HEADER  TITLE  FOR  Z(n)  FILE’ 

WRITE(**) 

READ(*:^00)  TITZ 
WRITE(*,*) 

WRITE(*,*)  ’SAVING  Z(t)  TO  OUTPUT  DATA  FILE’ 

WRrrE(3.100)  TITZ 
WRrrE(3.120)  N2 
WRITE(3,130)  TZ1 
WRITE(3,130)  T72 
DO  80  1=1, N2 
WRITE(3,130)  Z(l) 

80  CONTINUE 
CLOSE(3) 

C  QUERY  TO  DO  ANOTHER  RUN 
0 

95  CONTINUE 
WRITE(*,*) 

WRITE(MOO)  BELL 

WRITER, *)  ’DO  YOU  WISH  TO  CONVOLVE  ADDITIONAL  WAVEFORMS?’ 
WRITE(*,*) 

READ(MOO)  ANS 

IF((ANS.EQ.’Y’);OR.(ANS.EQ.’y’))THEN 
GO  TO  10 

ELSEIF((ANS.EQ.’N’).OR.(ANS.EQ.’n’))THEN 

GO  TO  99 

ELSE 

WRITE(^*)  ’TRY  AGAIN!  'Y'  OR  ’N'.’ 

GO  TO  95 
ENDIF 

99  CONTINUE 
WRITE(MOO)  BELL 

WRITE(*,*)  ’CONVOLUTION  PROGRAM  IS  OVER.  CHECK  NEW  FILE(S)’ 

100  FORMAT(A) 

120  F0RMAT(I5) 

130  FORMAT(E12.6) 

STOP 


APEENDIX  H.  INVERSE  ARMA  PROGRAM 


A.  PROGRAM  DESCRIPTION 

This  program  implements  an  inverse  ARMA  aspect  dependent  annihilation 
filter.  As  such,  it  is  an  adaptation  of  the  ARMA  signal  generator  in  Appendix  B.  It 
requires  input  of  the  and  6^.  coefficients  for  the  target  transfer  function  it  is 
intended  to  annihilate.  It  then  creates  the  inverse  impulse  response  and  convolves 
it  with  an  input  data  record.  Finally,  it  computes  the  average  energy  over  a  user 
specified  time  period  in  the  output. 

This  program  was  written  by  Prof.  M.  A.  Morgan  and  modified  by  Capt.  P.-C. 
Reddy. 


B.  PROGRAM  LISTING 

Program  lARMA 
C 

C  Reading  N-th  order  ARMA  filter  with  coefficients: 

C 

C  HR  a(l<)fork=1,N 

C  FIR  b(k)fork=0,L 

C 

C  Processing  y[n]  using  the  inverse  ARMA  filter  with  output 
C  z[n].  Inverse  ARMA  has  coefficients: 

C 

C  HR  ai(k)=-b{k)/b(0)  fork=1,L 

C  FIR  bi(k)=-a(k)/b(0)fork=1,N 

C  bi(0)=  1.0/b(0) 

C 

C  by  M.A.  Morgan  5/24/90 

C  Modified  by  P.C.  Reddy  5/27/90 

C 

REAL*8  y(0:2047), 2(0:2047) 

REAL*8  a(1 :30),b(0:32),ai(1 :30),bi(0:30) 
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CHARACTER*64  TITLE, TITL 
CHARACTER*T6  FNAME 
C  Entering;D.E.  Coefficients 

WRITE(*,*)  ’Enter  Filename  For  Recursive  a(k):’ 
WRITE(*,*)  ’Use  0  for  FIR  Filter’ 

READ(*,100)  FNAME 
IF(FNAME.EQ;’0’)  GOTO  15 
OPEN(2,FILE=FNAME) 

0  Assuming  that  a(k)  are  reversed  polarity  in  file. 
READ(2,110)N 
write(*,*)  ’Write  k,  a(k);’ 

DO  11  k=1,N 
READ(i*)  a(k) 
a(k)=-a(k) 

0  write(*,*)  k,a(k) 

11  CONTINUE 
CLOSE(2) 

15  CONTINUE 

WRITE(*,*)  ’Enter  Filename  for  Non-Recursive  b(k):’ 

READ(*.100)  FNAME 

OPEN(2,FILE=FNAME) 

READ(2,110)  L 

C  write(*,*)  ’Write  k,  b(k):’ 

DO  22  k=0,L 
READ(2,*)  b(k) 

C  write(*.*)  k,b(k) 

22  CONTINUE 
CLOSE(2) 

WRiTE(*,*)  ’Enter  Filename  for  y[n]  Plot  File;’ 

READ(*,100)  FNAME 

OPEN  (2, FILE= FNAME 

READ(2,100)  TITL 

READ  (2, 110)  NY 

READ(2,120)  XQ 

READ(2,T20)  XQ 

DO  24  m=0,NY-1 

READ(2,120)  y(m) 

24  CONTINUE 
CLOSE(2) 

WRITE(*,*) 

C  Setting  up  Plot  File  for  z[n] 

WRITE(*,*)  ’Enter  Number  of  Points  for  z[n):  (,le,  2048)’ 
READ(*.*)  NZ 

WRITE(*,*)  ’Enter  time  interval  (nanoseconds)’ 

READ(*,*)  Tmax 

Tmin=0.0 

WRITE(*,*)  ’Enter  filename  for  z[n]  response  output’ 
READ(*,100)  FNAME 

TITLE=  ’Inverse  ARMA  Filter  z[n]  Response’ 
OPEN(1,FILE=FNAME) 
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WRITE(1,100)=TITLE 
WRITE(1,110)  NZ 
WRITE(1,120)fTmin 
VyRITE(il20)Tmax 

C  Generating  lnverse;Filter  Coefficients 

C  Assuming  b(0)  .ne.  O  !!! 

Li=N 

Ni=L 

DO  50  k=1,Ni 
50  ai(k)=-b(k)/b(0) 

bi(0)=1.0D00/b(0) 

DO  55  k=1,U: 

55  bi(k)=-a(k)/b(0) 

WRITE(*,*) 

C  Initializing  then  iterating  to  form  z[m] 

z(b)=bi(0)*y(0) 
write(1,120)  2(0) 

DO  77  m=1,NZ-1 

z(m)=0.0 

Lmax=MIN(m,Li) 

DO  60  k=0,Lmax 
IF((m-k):GE.NY)  GO  TO  60 
z(m)=z(m)+bi{k)*y(m-k) 

60  Continue 

Nmax=MIN(m,Ni) 

DO  66  k=1,Nmax 
66  z(m)=z(m)+ai(k)*2(m-k) 

C  write(* ,*)  z(m) 

WR1TE(1,120)  z(m) 

77  CONTINUE 
CLOSE(1) 

C 

C  CALCULATING  LATE-TIME  ENERGY 


WRITE(* ,*)  ’ENTER  TL  (Nsec)  FOR  LATE-TIME  ENERGY  CALCUUTION' 
WRITE(**) 

READ(* ,*)  TLATE 

WRITE(* ,*)  ’ENTER  TE  (Nsec)  FOR  LATE-TIME  ENERGY  CALCUUTION’ 
WRITE(*,*) 

READ(* ,*)  TE 

NUTE=NINT((rLATE/Tmax)*NZ) + 1 
NE=NINT((TE^max)*NZ)+1 
ENERGY=0.0 
DO  90  N=NLATE,NE 
90  ENERGY=ENERGY+Z(N)*Z(N) 

ENERGY=ENERGY/(NE-NLATE+1) 

WRITE(*,*) 

WRITE(*,*)  ’LATE  TIME  ENERGY  =’,  ENERGY 
WRITE(* ,*)  Thats  All  Folks  !’ 
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100  FORMAT(A) 

110  F0RMAT(I5) 

120  FORM  At(E1 2.6) 
STOP 
END 
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